********> bugfix.16 Author: Ben Roberts and Gustavo Seabra Quantum Theory Project - University of Florida Date: Feb 11, 2009 Programs: sander.PUPIL Description: This bugfix only influences the sander.PUPIL interface. A number of small fixes to the sander/PUPIL interface are included here: - Some spelling corrections - Reduced the amount of information printed - Improved the behavior when exiting sander - Neighbor list is now updated every step if using PUPIL, to account for possible changes in the QM zone Fix: to apply this bugfix, save this file in $AMBERHOME and do the following: $ cd $AMBERHOME $ patch -p0 -N -r patch_rejects < bugfix.16 =========================================================== --- src/sander/force.f 2009-02-10 16:52:42.000000000 -0500 +++ src/sander/force.f 2009-02-11 12:08:29.000000000 -0500 @@ -219,6 +219,10 @@ if( ipimd>0.or.ineb>0) nrg_all(1:nbead)=0.d0 #endif + +! This shuld only happen if *not* using PUPIL, +! because PUPIL resets the nb list itself +#ifndef PUPIL_SUPPORT if( igb == 0 .and. iyammp == 0 ) then ! (for GB: do all nonbondeds together below) @@ -234,7 +238,8 @@ call timer_stop(TIME_LIST) call timer_stop(TIME_NONBON) end if - +#endif + #ifndef DISABLE_NCSU # ifdef MPI call ncsu_remember_initremd(rem.gt.0.and.mdloop.eq.0) @@ -359,7 +364,7 @@ ! Getting the Quantum forces with PUPIL package !***************************************************** -! Actualizing the simulation cell if there is any change +! Reconstructing the simulation cell if there is any change ! call inipupcell(natms,qcell,cell,xxx,yyy,zzz) do iPup=1,3 !vector loop do jPup=1,3 !Component loop @@ -371,7 +376,7 @@ qcell(11) = 0.0d0 qcell(12) = 0.0d0 -! temporary vector to wrap the real coordinates to pass throught +! temporary vector to wrap the real coordinates to pass through ! PUPIL interface. call get_stack(l_puptmp,3*natom,routine) if(.not. rstack_ok)then @@ -389,7 +394,7 @@ if(ifbox == 2) call wrap_to(nspm,ix(i70),r_stack(l_puptmp),box) end if -! write(6,*) ' Starting to update data PUPIL structures.' +! write(6,*) ' Updating PUPIL data structures.' ! Preparing the coordinates,velocity and classic forces ! to get quantum force @@ -488,13 +493,13 @@ enddo endif -! Getting the quantum forces for an specific quantumn domain +! Getting the quantum forces for a specific quantumn domain pupStep = pupStep + 1 puperror = 0 pupLevelData = 3 call getquantumforces(natom,pupLevelData,pupStep,puperror,qcdata,qcell) if (puperror.ne.0) then - write (6,*) "Error getting quantum forces." + write (6,*) "Fatal Error: Error getting quantum forces." call mexit(6,1) endif ! Quantum energy treatment.... @@ -565,11 +570,29 @@ enddo call deleting_qm_atoms() - + qsetup = .true. ! Setting as current quantum zone pupQZchange = 0 endif + ! For PUPIL, rebuild the neighbour list + ! and zero the charges on QM atoms at every step, + ! because the QM atoms list may have changed + if ( igb == 0 .and. iyammp == 0 ) then + + ! (for GB: do all nonbondeds together below) + call timer_start(TIME_NONBON) + call timer_start(TIME_LIST) + !do_list_update=.true. + call nonbond_list(x,ix(i04),ix(i06),ix(i08),ix(i10), & + ntypes,natom/am_nbead,xx,ix,ipairs,ntnb, & + ix(ibellygp),belly,newbalance,cn1, & + xx(lvel),xx(lvel2),ntp,xx(l45), qsetup, & + do_list_update) + !call qm_zero_charges(x(L15)) + call timer_stop(TIME_LIST) + call timer_stop(TIME_NONBON) + end if #endif /*PUPIL_SUPPORT*/ --- src/sander/mexit.f 2007-10-08 23:32:41.000000000 -0400 +++ src/sander/mexit.f 2009-02-10 16:51:14.000000000 -0500 @@ -36,9 +36,11 @@ #ifdef PUPIL_SUPPORT !jtc ========================= PUPIL INTERFACE ========================= ! Terminating the PUPIL corba Interface - puperror = 0 - call killcorbaintfc(puperror) - if(puperror /= 0) write(6,*) 'Error ending PUPIL CORBA Interface.' + if(pupactive) then + puperror = 0 + call killcorbaintfc(puperror) + if(puperror /= 0) write(6,*) 'Error ending PUPIL CORBA Interface.' + end if !jtc ========================= PUPIL INTERFACE ========================= #endif --- src/sander/pupildata.f 2007-12-10 11:45:23.000000000 -0500 +++ src/sander/pupildata.f 2009-02-10 16:51:14.000000000 -0500 @@ -1,6 +1,7 @@ #include "dprec.h" module pupildata + logical pupactive integer puperror,iPup,jPup,bs1,bs2,iresPup,l_puptmp integer pupLevelData ! Data to pass throught PUPIL Interface integer pupStep @@ -122,16 +123,16 @@ subroutine deleting_qm_atoms() use pupildata, x=>realStack , ix=>ixStack, ih=>ihStack - use parms, only: req,fmn + use parms, only: req, fmn #include "memory.h" #include "extra_pts.h" - + integer replicates integer i replicates = 1 -! Initializing the list of structrues for a new QM zone +! Initializing the list of structures for a new QM zone call init_extra_pts( & ix(iibh),ix(ijbh),ix(iicbh), & ix(iiba),ix(ijba),ix(iicba), & @@ -160,18 +161,18 @@ if(nphia.gt.0) call setdih(nphia,ix(i50),ix(i52),ix(i54),ix(i56), & ix(i58), ix(ibellygp)) ! remove dihedrals between QM atoms from list - if(pupQZChange .ne. 0) then - write(6,*) 'NUMBER OF H BONDS = ',nbonh - write(6,*) 'NUMBER OF BONDS = ',nbona - write(6,*) 'NUMBER OF H ANGLES = ',ntheth - write(6,*) 'NUMBER OF ANGLES = ',ntheta - write(6,*) 'NUMBER OF H DIHEDRALS = ',nphih - write(6,*) 'NUMBER OF DIHEDRALS = ',nphia - write(6,*) 'NUMBER OF NONBONDED 14 PAIRS = ',numnb14 - do i=1,numnb14 - write(6,*) ix(inb_14+(i-1)*2),ix(inb_14+(i-1)*2+1) - enddo - endif + !if(pupQZChange .ne. 0) then + ! write(6,*) 'NUMBER OF H BONDS = ',nbonh + ! write(6,*) 'NUMBER OF BONDS = ',nbona + ! write(6,*) 'NUMBER OF H ANGLES = ',ntheth + ! write(6,*) 'NUMBER OF ANGLES = ',ntheta + ! write(6,*) 'NUMBER OF H DIHEDRALS = ',nphih + ! write(6,*) 'NUMBER OF DIHEDRALS = ',nphia + ! write(6,*) 'NUMBER OF NONBONDED 14 PAIRS = ',numnb14 + ! do i=1,numnb14 + ! write(6,*) ix(inb_14+(i-1)*2),ix(inb_14+(i-1)*2+1) + ! enddo + !endif ! write(6,*) 'nbonh',nbonh,'iibh',iibh,'nbona',nbona,'ntheth',ntheth ! write(6,*) 'ntheta',ntheta,'nphih',nphih,'nphia',nphia --- src/sander/putvalues.f 2007-01-18 17:09:50.000000000 -0500 +++ src/sander/putvalues.f 2009-02-10 16:51:14.000000000 -0500 @@ -5,11 +5,11 @@ !C C !C PUPIL program ITR-MEDIUM PROJECT C !C C -!C Author: J. Torras C -!C Version: 0.0.8.9 C -!C Date: 09-08-2005 C -!C Proposal Date: 08-08-2005 C -!C Revising Date: 02-06-2006 C +!C Author: J. Torras C +!C Version: 0.0.8.9 C +!C Date: 09-08-2005 C +!C Proposal Date: 08-08-2005 C +!C Revising Date: 02-06-2006 C !C C !C torras@qtp.ufl.edu C !C C @@ -31,9 +31,10 @@ !C EXTERNAL CALCULATION OF FORCE C !C C !C NT TOTAL NUMBER OF ATOMS ASSOC. WITH EXTERNAL CALCULATION OF FORCES C -!C QZ IF QZ IS DIFFERENT OF ZERO, THE QUANTUM ZONE HAS BEEN CHANGED C -!C THE NEW DEFINITIONS OF QUANTUM ZONE IS GIVEN IN A C -!C A INTEGER VECTOR OF DIMENSION N*2 WITH THE ATOM NUMBER AND ATOM ROLE C +!C QZ IF QZ IS NON-ZERO, THE QUANTUM ZONE HAS BEEN CHANGED. C +!C THE NEW DEFINITION OF THE QUANTUM ZONE IS GIVEN IN AN... C +!C A ...INTEGER VECTOR OF DIMENSION 2*NT WITH THE ATOM NUMBER AND ATOM C +!C ROLE C !C ATOM_NUM1,ATOM_ROLE1, ... , ATOM_NUMn,ATOM_ROLEn C !C Atom Role: 0 -> Quantum atom C !C 1 -> Pseudoatom C @@ -54,6 +55,7 @@ qmmm_scratch #include "memory.h" +#include "debug.h" INTEGER a (nt*2),nt,qz _REAL_ f (nt*3),e @@ -64,13 +66,18 @@ ! Initial considerations under QZ Change...... pupQZchange = qz + +! If we are debugging, assume the QM zone has changed. + if ( do_debugf /= 0 ) pupQZchange = 1 + if (pupQZchange .ne.0) then ! QM-zone changed - write(6,*) ' NEW QUANTUM ZONE ==> ',nt,' PARTICLES' + if (qmmm_nml%verbosity > 2) & + write(6,"(a16,1x,i6,1x,a23)") 'Searching within', nt, 'particles for QM atoms.' ! Restoring charges to the previous set of quantum particles... do i=1,pupqatoms - write(6,"(a12,2x,i6,2x,a4,f12.6,2x,a4,f12.6)") & - 'Charge QUAT=',pupqlist(i),'Old=',realStack(L15+pupqlist(i)-1),& - 'New=',pupchg(pupqlist(i)) + ! write(6,"(a12,1x,i6,3x,a12,1x,f12.6,3x,a12,1x,f12.6)") & + ! 'Atom number:',pupqlist(i),'Old charge =',realStack(L15+pupqlist(i)-1),& + ! 'New charge =',pupchg(pupqlist(i)) realStack(L15+pupqlist(i)-1) = pupchg(pupqlist(i)) enddo @@ -78,14 +85,15 @@ do i=1,natom pupmask(i) = 0 enddo - + ! Counting and getting real quantum particles totQM = 0 do i=1,nt pupqlist(i) = abs(a((i-1)*2+1)) if(a((i-1)*2+2).eq.0) then pupmask(a((i-1)*2+1)) = 1 - write(6,*) 'NEW QA', a((i-1)*2+1) + if (qmmm_nml%verbosity > 2) & + write(6,"(a7,i6,a1,1x,i6)") 'QM atom', a((i-1)*2+1), ':', i totQM = totQM + 1 endif enddo @@ -105,16 +113,16 @@ REQUIRE(ier == 0) endif -! Initializing qmm_nml and qmmm_struct structures if defined.... +! Initializing qmm_nml and qmmm_struct structures if these are present. igb = 0 ier = 0 - qmmm_struct%nquant = totQM ! number of current quantum atoms - qmmm_struct%nlink = 0 ! no link atoms - qmmm_nml%qm_pme = .false. ! no calculate QM-QM long interactions ? - qm2_struct%calc_mchg_scf = .false. ! never calculate Mulliken charges - qmmm_nml%qmshake = 0 ! has no sense QM shaking - qmmm_nml%qm_ewald = 0 ! No QM-QM ewald interaction + qmmm_struct%nquant = totQM ! Current number of QM atoms + qmmm_struct%nlink = 0 ! Number of link atoms (default: 0) + qmmm_nml%qm_pme = .false. ! Use PME to compute long-range QM-QM electrostatics? (default: No) + qm2_struct%calc_mchg_scf = .false. ! Compute Mulliken charges? (default: No) + qmmm_nml%qmshake = 0 ! Use SHAKE? (default: No) + qmmm_nml%qm_ewald = 0 ! Use Ewald interpolation for long-range QM-QM electrostatics? (default: No) allocate ( qmmm_scratch%qm_real_scratch(3*natom),stat=ier) REQUIRE(ier == 0) allocate ( qmmm_struct%iqm_atomic_numbers(totQM),stat=ier) @@ -128,9 +136,11 @@ qmmm_struct%atom_mask(1:natom) = .false. qmmm_struct%qm_resp_charges(1:totQM) = zero - endif ! if (pupQZchange .ne.0) ! End of QUANTUM ZONE CHANGED - -! Assing the number of particles with quantum force contribution + endif + +! From here on, do whether or not the QM zone has changed + +! Determine the number of particles in the PUPIL QM region pupqatoms = nt totQM = 0 @@ -138,30 +148,40 @@ bs1 = (pupqlist(i)-1)*3 bs2 = (i-1)*3 do j=1,3 -! Assing directly the quantum force +! Calculate the quantum force qfpup(bs1+j) = f(bs2+j)*EV_TO_KCAL enddo ! write(6,"(a22,2x,i4,3(2x,e16.10),2x,a8,2x,f10.7)") 'Quantum Particle ==>F:',pupqlist(i),(f(bs2+k)*EV_TO_KCAL,k=1,3),'Charge =',chg(i) - -! Managing the quantum list - if(qmmm_nml%ifqnt .and. a((i-1)*2+2).eq.0 .and. (pupQZchange.ne.0) ) then + + ! Manage the quantum list + if (qmmm_nml%ifqnt .and. a((i-1)*2+2) .eq. 0) then + if (pupQZchange .ne. 0) then totQM = totQM + 1 qmmm_struct%iqm_atomic_numbers(totQM) = pupatm(pupqlist(i)) qmmm_nml%iqmatoms(totQM) = pupqlist(i) qmmm_struct%atom_mask(pupqlist(i)) = .true. -! Putting quantum atom charges to zero + ! Putting quantum atom charges to zero + ! BPR - try using Mulliken charges instead realStack(L15+pupqlist(i)-1) = zero -! write(6,"(a12,2x,i6,2x,a4,f12.6,2x,a4,f12.6)") 'Charge Atom=',pupqlist(i),'Old=',realStack(L15+pupqlist(i)-1),'New=',chg(i)*AMBER_ELECTROSTATIC -! realStack(L15+pupqlist(i)-1) = chg(i) * AMBER_ELECTROSTATIC - endif - + endif + ! This instruction was originally inside the pupQZchange .ne. 0 loop. + ! But, of course, Mulliken charges can change even if the atoms in the QM region do not. + !write(6,*) 'Updating charges on QM atoms' + !write(6,*) '----------------------------' + !write(6,"(a12,2x,i6,2x,a4,f12.6,2x,a4,f12.6)") 'Charge Atom=',pupqlist(i),'Old=',realStack(L15+pupqlist(i)-1),'QM Charge=',chg(i)*AMBER_ELECTROSTATIC + !realStack(L15+pupqlist(i)-1) = chg(i) * AMBER_ELECTROSTATIC + endif enddo -! QM Energy - qmEnergy = E * EV_TO_KCAL - write(6,"(a11,e16.10)") 'QM ENERGY =',qmEnergy - RETURN - END + ! Print the total QM (SCF) energy to mdout + qmEnergy = E * EV_TO_KCAL + if (qmmm_nml%verbosity > 2) then + write(6,"(a12,e18.10)") 'QM ENERGY = ', qmEnergy + write(6,"") + write(6,"") + endif + RETURN +END !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC --- src/sander/sander.f 2009-02-10 16:52:42.000000000 -0500 +++ src/sander/sander.f 2009-02-11 12:18:35.000000000 -0500 @@ -179,6 +179,7 @@ write(6,*) 'Error in the PUPIL interface initialization.' call mexit(6,1) endif + pupactive = .true. write(6,*) 'PUPIL CORBA Interface initialized.' #endif !jtc ========================= PUPIL INTERFACE ========================= @@ -373,7 +374,7 @@ ! Getting all the initial charges.... pupchg(iPup) = x(L15+iPup-1) - write(6,*) 'Atom num.',iPup,' Label ', keyMM(iPup), 'Charge', pupchg(iPup) + ! write(6,*) 'Atom num.',iPup,' Label ', keyMM(iPup), 'Charge', pupchg(iPup) do jPup=1,3 qfpup(bs1+jPup) = 0.0d0 @@ -396,9 +397,9 @@ ! Submitting the Residue Pointer vector to PUPIL write(6,*) 'Number of residues = ',nres,' Number of atoms= ',natom - do iPup=1,nres - write(6,*) 'Residue ',iPup,keyres(iPup),pupres(iPup) - enddo + !do iPup=1,nres + ! write(6,*) 'Residue ',iPup,keyres(iPup),pupres(iPup) + !enddo puperror = 0 call putresiduetypes(nres,puperror,pupres,keyres) if (puperror .ne. 0) then @@ -1139,6 +1140,7 @@ if(puperror /= 0) then write(6,*) 'Error ending PUPIL CORBA Interface.' endif + pupactive = .false. write(6, '(a)') 'PUPIL CORBA Interface finalized.' #endif !! jtc ========================= PUPIL INTERFACE =========================