C C $Id: edriver.F,v 1.15 1999/03/02 20:02:25 arjan Exp arjan $ C C------------------------------------------------------------------------ subroutine edriver(iflag, iflaggr, nat, $ xyzroar, gradroar, ierror,gradpme, symbol2) use se_corrections_module, only : Apply_se_corrections implicit double precision (a-h,o-z) #include "divcon.dim" #include "divcon.h" #include "constants.h" #ifdef MPI #include "mpif.h" #endif dimension xyzroar(*), gradroar(*),gradpme(*) dimension gradpme1(maxatm), gradpme2(maxatm),grad1(3,maxatm) dimension grade(3*maxatm) C locals: logical newsub integer ierror integer tmp_mpi_data ! Antoine MARION :: 2013-12-10 Start dimension gradSEcorrections(3,3*maxatm) dimension vir_SE(3) character*4 symbol2(nat) double precision E_SEcorrections, E_SEpep logical SE_corrections double precision vir_se_corrections SE_corrections = .true. ! SE_corrections = .false. ! Antoine MARION :: 2013-12-10 End CCC-----------------------------------------------------------------CCC C this part is not parallelized, run it by the master node ierror = 0 if (iflag.eq.0) then C . set the convergence criteria i = index(keywrd,'ECRIT=') if (i.ne.0) then call rdnum(keywrd,i,eecrt,ierror) endif i = index(keywrd,'DCRIT=') if (i.ne.0) then call rdnum(keywrd,i,dencrt,ierror) endif endif if (iflag.eq.1) then newsub = .true. else newsub = .false. endif C print *, 'newsub ', newsub, iflag if (iflag.eq.0) then IF (rotat.and.rotat2) THEN continue ELSE call setup(ierror) if (SE_corrections) then call setupSEcorrections(symbol2,ierror) endif ENDIF if(ierror.ne.0) go to 1000 endif C if (boxroar(1).gt.0.1) then C dbox(1) = boxroar(1) C dbox(2) = boxroar(2) C dbox(3) = boxroar(3) C endif C dhalf(1) = 0.5d0*dbox(1) C dhalf(2) = 0.5d0*dbox(2) C dhalf(3) = 0.5d0*dbox(3) if (nat.ne.natoms) then write(*,'(" ERROR: THE NUMBER OF ATOMS FROM ROAR AND DIVCON", & " IS NOT IDENTICAL")') ierror = 1 return endif k3 = 3 k2 = 2 k1 = 1 do i=1,natoms xyz(1,i) = xyzroar(k1) xyz(2,i) = xyzroar(k2) xyz(3,i) = xyzroar(k3) k1 = k1 + 3 k2 = k2 + 3 k3 = k3 + 3 enddo if (iflag.eq.0) then if(index(keywrd,'OPT=').ne.0)then write(iout,'(/" ERROR: NO GEOMETRY OPTIMIZATION ", & "WITHIN A MD SIMULATION")') ierror =1 return wrtedmx = .false. if (pme) then write(iout,'(/" ERROR: GEOMETRY OPTIMIZATION WITH ", & "PME IS NOT IMPLEMENTED")') ierror = 1 goto 1000 elseif (inter) then write(iout,'(/" ERROR: INTER KEYWORD IS NOT ALLOWED ", & "FOR GEOMETRY OPTIMIZATION")') ierror = 1 goto 1000 endif if(index(keywrd,'TESTRUN').ne.0) return call geoopt(eheat,ierror) c call atmchg else if(ierror.ne.0) goto 1000 if (index(keywrd,"PUSH").ne.0) then write(iout,'(/" ERROR: NO PUSH WITHIN A MD SIMULATION")') ierror = 1 return call push endif C . setup for PME endif endif nscf = 0 call setbox(ierror) if (iflag.eq.1) then if (pme) then call pme_setup call pme_calcb call pme_calctheta C-RDC write(iout,'( C-RDC $ /" A CLASSICAL SMOOTH PARTICLE MESH EWALD ", C-RDC & "SUM (PME) IS USED TO", C-RDC & /" CALCULATE THE LONG RANGE ", C-RDC $ "ELECTROSTATIC INTERACTION.", C-RDC & //" PME-PARAMETERS:",/" --------------")') C-RDC write(iout,'(" NUMBER OF GRIDS IN THE X-DIRECTION: ", C-RDC & I4,/" NUMBER OF GRIDS IN THE Y-DIRECTION: ", I4, C-RDC & /" NUMBER OF GRIDS IN THE Z-DIRECTION: ", I4, C-RDC & /" ORDER OF CARDINAL B-SPLINE INTERPOLATION: ", C-RDC $ I4, C-RDC & /" EXPONENTIAL BETA-FACTOR: ", F10.6)') C-RDC & k1pme,k2pme,k3pme,nspline,betapme ! if (cm1) then C-RDC write(iout,'(/" PME WILL USE CM1 CHARGES")') ! elseif (cm2) then C-RDC write(iout,'(/" PME WILL USE CM2 CHARGES")') ! else C-RDC write(iout,'(/" PME WILL USE MULLIKEN CHARGES")') ! endif ! if (inter.and.recipintr) then C-RDC write(iout,'(/" ALL INTRAMOLECULAR PME CONTRIBUTIONS", C-RDC & " TO THE GRADIENT ARE NEGLECTED,", C-RDC & /" EXCEPT FOR THE INTRAMOLECULAR CONTRIBUTIONS", C-RDC & " OF THE RECIPROCAL ENERGY")') ! elseif (inter) then C-RDC write(iout,'(/" ALL INTRAMOLECULAR PME CONTRIBUTIONS", C-RDC & " TO THE GRADIENT ARE NEGLECTED")') ! else C-RDC write(iout,'(/" BOTH INTRAMOLECULAR AND INTERMOLECULAR", C-RDC & " PME CONTRIBUTIONS", C-RDC & /" ARE ADDED TO THE GRADIENT")') ! endif endif endif if (iflag.eq.0) return if (newsub) then CALL GENSUB(IERROR) IF(IERROR.NE.0) RETURN endif if(index(keywrd,'TESTRUN').ne.0) return call energy(newsub,eheat,ierror) ! Antoine MARION :: 2013-12-10 Start if (SE_corrections) then call Apply_se_corrections(XYZ,E_SEcorretions,E_SEpep, & gradSEcorrections,vir_SE) ecore = ecore + E_SEcorretions etot = etot + E_SEcorretions eheat = eheat + (E_SEcorretions + E_SEpep)*ev2kcal endif ! Antoine MARION :: 2013-12-10 End if(ierror.ne.0) goto 1000 call atmchg if (pme) then if (gradient) then c if (inter) then c call gcartmc(grad) c else call gcart(grad1) ii=-3 do i=1,natoms ii=ii+3 grade(ii+1)=grad1(1,i) grade(ii+2)=grad1(2,i) grade(ii+3)=grad1(3,i) enddo c endif endif do i=1,3*natoms gradpme1(i)=0.0d0 gradpme2(i)=0.0d0 enddo c if (inter) then c call pme_directmc(edir,eself,ecoul,grad) c else call pme_direct2(edir,eself,ecoul,gradpme1) c endif c if (recipintr) then call pme_recip(erec,gradpme2) c else c call pme_recipmc(erec,grad) c endif do i=1,3*natoms gradpme(i)=gradpme1(i)+gradpme2(i) grade(i)=grade(i)+gradpme(i) enddo C . add PME energies to total energy c elr = edir + eself + erec - ecoul elr = edir + erec etot = etot + elr C . add PME contribution to heat of formation ! eheat = eheat + 23.061*elr eheat = eheat + eV2kcal*elr do i=1,3 vir(i)=vir(i)+eV2kcal*(-elr-ecoul)/3.0d0 enddo c do i=1,3*natoms c gradpme(i)=23.061*gradpme(i) c enddo endif if (scrf) then C-RDC write(*,'(" ERROR: SCRF NOT WITHIN A MD SIMULATION")') ierror = 1 return C-RDC write(iout,'(//" FINAL QUANTITIES ", C-RDC $ "(SOLVENT NOT INCLUDED):", C-RDC . /" ---------------------------------------")') else C-RDC write(iout,'(//" FINAL QUANTITIES:", C-RDC . /" -----------------")') endif C-RDC write(iout,'(/" NUMBER OF FULLY DIAGONALIZED SUBSYSTEMS: ",I6, C-RDC & /" NUMBER OF FROZEN SUBSYSTEMS: ",I6, C-RDC & /" NUMBER OF SCF CALCULATIONS = ",i6)') C-RDC & nfull, nfroz, nscf C-RDC write(iout,'(/" ELECTRONIC ENERGY = ",f25.8," EV", C-RDC . /" CORE-CORE REPULSIONS = ",f25.8," EV")') C-RDC . eelect,ecore C-RDC if (pme) write(iout,'(" PME DIRECT ENERGY = ",f25.8, C-RDC & " EV",/" PME RECIPROCAL ENERGY = ",f25.8," EV", C-RDC & /" PME SELF ENERGY = ",f25.8," EV", C-RDC & /" CLAS. COULOMB ENERGY: = ",f25.8," eV", C-RDC & /" TOTAL LONGE RANGE = ",f25.8," eV")') C-RDC & edir,erec, eself, ecoul, elr C-RDC write(iout,'(" TOTAL ENERGY = ",f25.8," EV", C-RDC . /" HEAT OF FORMATION = ",f25.8," KCAL/MOL")') C-RDC . etot,eheat if (setch) then do i=1,ncores C-RDC write(*,'(" FERMI ENERGY OF CLUSTER GROUP ",I4," = ", C-RDC & f25.8," EV")') C-RDC & i, efermi(i) enddo else C-RDC write(*,'(" FERMI ENERGY = ",f25.8," EV")') C-RDC & efermi(1) endif C-RDC IF(ionpot) write (iout, C-RDC $ '(" IONIZATION POTENTIAL = ",f25.8," EV")') C-RDC $ -eval(nelec/2) call homolumo if(screen) then if (scrf) then write(iout,'(" DIVCON FINAL QUANTITIES ", $ "(SOLVENT NOT INCLUDED):", . /" -----------------")') else write(iout,'(" DIVCON FINAL QUANTITIES:", . /" -----------------")') endif #ifdef MPI call mpi_allreduce(nfull,tmp_mpi_data,1,MPI_INTEGER,MPI_SUM, . commsebomd,ier) nfull = tmp_mpi_data #endif write(iout, $ '(" NUMBER OF FULLY DIAGONALIZED SUBSYSTEMS: ",I6, & /" NUMBER OF FROZEN SUBSYSTEMS: ",I6, & /" NUMBER OF SCF CALCULATIONS = ",i6, & /" NUMBER OF DIAGONALIZATIONS = " ,i6)') & nfull, nfroz, nscf, ndiag write(iout,'(/" ELECTRONIC ENERGY = ",f25.8," EV", . /" CORE-CORE REPULSIONS = ",f25.8," EV")') . eelect,ecore if (pme) write(iout,'(" PME DIRECT ENERGY = ",f25.8, & " EV",/" PME RECIPROCAL ENERGY = ",f25.8," EV", & /" PME SELF ENERGY = ",f25.8," EV", & /" CLAS. COULOMB ENERGY = ",f25.8," EV", & /" TOTAL LONG RANGE = ",f25.8," EV")') & edir,erec,eself,ecoul, elr ! & pmedir,pmerec,pmeself,ecoul, elr write(iout,'(" TOTAL ENERGY = ",f25.8," EV", . /" HEAT OF FORMATION = ",f25.8," KCAL/MOL", . /" VIRIAL = ",f25.8)') . etot,eheat,vir(4) if (setch) then do i=1,ncores write(iout,'(" FERMI ENERGY OF CLUSTER GROUP ",I4," = ", & f25.8," EV")') & i, efermi(i) enddo else write(iout,'(" FERMI ENERGY = ",f25.8," EV")') & efermi(1) endif IF(ionpot) write (iout, $ '(" IONIZATION POTENTIAL = ",f25.8," EV")') $ -eval(nelec/2) write(iout,'(" -----------------")') endif IF(SCRF) THEN if (scrfin) then C-RDC write(iout,'(/" FILLING THE CAVITIES WITH ", C-RDC $ "A DIELECTRIC CONTINUUM"/)') endif time_SCF = 0.D0 time_PB = 0.D0 if(.NOT.SOLVCALC) then C . No initial surface charges were read in or restart job in PB solver SOLVCALC = .TRUE. ! Turn solvation calculation ON E_VACC = EHEAT endif C . Backup vaccuum density matrix call system("mv divcon.dmx divcon.dmx.vac") C-RDC if (screen) C-RDC $ write(iscr, C-RDC $ '(" ... backup the vaccuum density matrix file")') C C . CALCULATE SCRF ENERGY C call ENERGY(newsub,EHEAT,IERROR) if(ierror.NE.0) return C-RDC write(iout,'(//" SOLUTION QUANTITIES:", C-RDC . /" -------------------")') C-RDC write(iout, C-RDC $ '(/" NUMBER OF FULLY DIAGONALIZED SUBSYSTEMS: ",I6, C-RDC & /" NUMBER OF FROZEN SUBSYSTEMS: ",I6, C-RDC & /" NUMBER OF SCF CALCULATIONS = ",i6)') C-RDC & nfull, nfroz, nscf C-RDC write(iout, C-RDC . '(/" ELECTRONIC ENERGY = ",f25.8," EV", C-RDC . /" CORE-CORE REPULSIONS = ",f25.8," EV")') C-RDC . eelect,ecore C-RDC if (pme) write(iout,'(" PME DIRECT ENERGY = ",f25.8, C-RDC & " EV",/" PME RECIPROCAL ENERGY = ",f25.8," EV", C-RDC & /" PME SELF ENERGY = ",f25.8," EV", C-RDC & /" CLAS. COULOMB ENERGY: = ",f25.8," eV", C-RDC & /" TOTAL LONGE RANGE = ",f25.8," eV")') C-RDC & edir,erec, eself, ecoul, elr C-RDC write(iout,'(" TOTAL ENERGY = ",f25.8," EV", C-RDC . /" HEAT OF FORMATION = ",f25.8," KCAL/MOL")') C-RDC . etot,eheat if (setch) then do i=1,ncores C-RDC write(iout, C-RDC $ '(" FERMI ENERGY OF CLUSTER GROUP ",I4," = ", C-RDC & f25.8," EV")') C-RDC & i, efermi(i) enddo else C-RDC write(iout,'(" FERMI ENERGY = ",f25.8," EV")') C-RDC & efermi(1) endif C-RDC if (ionpot) write (iout, C-RDC $ '(" IONIZATION POTENTIAL = ",f25.8," EV")') C-RDC $ -eval(nelec/2) call homolumo if (screen) then C-RDC write(iscr,'(//" SOLUTION QUANTITIES:", C-RDC . /" -------------------")') C-RDC write(iscr, C-RDC $ '(/" NUMBER OF FULLY DIAGONALIZED SUBSYSTEMS: ",I6, C-RDC & /" NUMBER OF FROZEN SUBSYSTEMS: ",I6, C-RDC & /" NUMBER OF SCF CALCULATIONS = ",i6)') C-RDC & nfull, nfroz, nscf C-RDC write(iscr, C-RDC . '(/" ELECTRONIC ENERGY = ",f25.8," EV", C-RDC . /" CORE-CORE REPULSIONS = ",f25.8," EV")') C-RDC . eelect,ecore C-RDC if (pme) write(iscr,'(" PME DIRECT ENERGY = ",f25.8, C-RDC & " EV",/" PME RECIPROCAL ENERGY = ",f25.8," EV", C-RDC & /" PME SELF ENERGY = ",f25.8," EV", C-RDC & /" CLAS. COULOMB ENERGY: = ",f25.8," eV", C-RDC & /" TOTAL LONGE RANGE = ",f25.8," eV")') C-RDC & edir,erec, eself, ecoul, elr C-RDC write(iscr,'(" TOTAL ENERGY = ",f25.8," EV", C-RDC . /" HEAT OF FORMATION = ",f25.8," KCAL/MOL")') C-RDC . etot,eheat if (setch) then do i=1,ncores C-RDC write(iscr, C-RDC $ '(" FERMI ENERGY OF CLUSTER GROUP ",I4," = ", C-RDC & f25.8," EV")') C-RDC & i, efermi(i) enddo else C-RDC write(iscr,'(" FERMI ENERGY = ",f25.8," EV")') C-RDC & efermi(1) endif C-RDC if (ionpot) write (iscr, C-RDC $ '(" IONIZATION POTENTIAL = ",f25.8," EV")') C-RDC $ -eval(nelec/2) endif C C . CALCULATE THE REORGANIZATION ENERGY with the Fock matrix C . without solvent contribution but with the same density matrix C . got in the last SCRF cycle. C . (change of the gas phase wave function due to solvation) C last_SCRF = .TRUE. SOLVCALC = .False. C-RDC if (screen) C-RDC $ WRITE(ISCR,'(/" PAST SCRF: ...CALCULATE TOTAL ENERGY", C-RDC . " WITHOUT SOLVENT CONTRIBUTION")') C-RDC WRITE(iout,'(/" PAST SCRF: ...CALCULATE TOTAL ENERGY", C-RDC . " WITHOUT SOLVENT CONTRIBUTION")') call ENERGY(newsub,EHEAT,IERROR) Gwfd = EHEAT- E_VACC time_tot = time_SCF + time_PB C-RDC write(iout,'(/" SUMMARY OF SCRF CALCULATION: ", C-RDC $ "(FREE ENERGIES [KCAL/MOL]): ", C-RDC . /" (1) REACTION FIELD (VACCUUM WAVE FUNCTION) = ", C-RDC $ f10.2," ", C-RDC . /" (2) REACTION FIELD (SOLUTION WAVE FUNCTION) = ", C-RDC $ f10.2," ", C-RDC . /" (3) WAVE FUNCTION DISTORTION = ", C-RDC $ f10.2," ", C-RDC . /" (4) POLARIZATION OF SOLUTE CHARGE = ", C-RDC . f10.2," (2-1+3)", C-RDC . /" (5) ELECTROSTATIC IN SOLUTION = ", C-RDC . f10.2," (2+3)", C-RDC . /" (6) CAVITY/DISPERSION = ", C-RDC $ f10.2," ", C-RDC . /" (7) SOLVATION = ", C-RDC . f10.2," (5+6)", C-RDC . /" SCF CPU (SEC) = ",f10.2," (", C-RDC . f5.1,"%)", C-RDC . /" PB CPU (SEC) = ",f10.2," (", C-RDC . f5.1,"%)")') C-RDC $ grf_vac, grf, gwfd, GRF-GRF_VAC+GWFD, grf+gwfd, C-RDC $ gnp, grf+gnp+gwfd, C-RDC $ time_scf,100*TIME_SCF/TIME_TOT, C-RDC $ time_pb,100*TIME_PB/TIME_TOT C-RDC if (screen) C-RDC $ write(iscr,'(/" SUMMARY OF SCRF CALCULATION: ", C-RDC $ "(FREE ENERGIES [KCAL/MOL]): ", C-RDC . /" (1) REACTION FIELD (VACCUUM WAVE FUNCTION) = ", C-RDC $ f10.2," ", C-RDC . /" (2) REACTION FIELD (SOLUTION WAVE FUNCTION) = ", C-RDC $ f10.2," ", C-RDC . /" (3) WAVE FUNCTION DISTORTION = ", C-RDC $ f10.2," ", C-RDC . /" (4) POLARIZATION OF SOLUTE CHARGE = ", C-RDC . f10.2," (2-1+3)", C-RDC . /" (5) ELECTROSTATIC IN SOLUTION = ", C-RDC . f10.2," (2+3)", C-RDC . /" (6) CAVITY/DISPERSION = ", C-RDC $ f10.2," ", C-RDC . /" (7) SOLVATION = ", C-RDC . f10.2," (5+6)", C-RDC . /" SCF CPU (SEC) = ",f10.2," (", C-RDC . f5.1,"%)", C-RDC . /" PB CPU (SEC) = ",f10.2," (", C-RDC . f5.1,"%)")') C-RDC $ grf_vac, grf, gwfd, GRF-GRF_VAC+GWFD, grf+gwfd, C-RDC $ gnp, grf+gnp+gwfd, C-RDC $ time_scf,100*TIME_SCF/TIME_TOT, C-RDC $ time_pb,100*TIME_PB/TIME_TOT ENDIF if(GRADIENT)then call setopt if (.not.pme) then c if (inter) then c call gcartmc(grad) c else call gcart(grad1) ii=-3 do i=1,natoms ii=ii+3 grade(ii+1)=grad1(1,i) grade(ii+2)=grad1(2,i) grade(ii+3)=grad1(3,i) enddo c do i=1,natoms c ii=ii+3 c grad(ii+1)=grad1(1,i) c grad(ii+2)=grad1(2,i) c grad(ii+3)=grad1(3,i) c enddo c endif endif ! Antoine MARION :: 2013-12-10 Start if (SE_corrections) then ii=-3 do i=1,natoms ii=ii+3 grade(ii+1)=grade(ii+1)+gradSEcorrections(1,i) grade(ii+2)=grade(ii+2)+gradSEcorrections(2,i) grade(ii+3)=grade(ii+3)+gradSEcorrections(3,i) enddo vir_se_corrections = vir_se(1)+vir_se(2)+vir_se(3) endif ! Antoine MARION :: 2013-12-10 End vir(4)=vir(1)+vir(2)+vir(3)+vir_se_corrections call dograd(grade,gnorm,ierror) if (ierror.ne.0) goto 1000 endif call wrttims ccccccccccccccccccccccccccccccccccccccccccccc c call dovir(vir, grad) ccccccccccccccccccccccccccccccccccccccccccccc C-RDC write(*,*) "edir eself ecoul erec ", edir, eself, ecoul, erec c if (iflaggr.eq.1) then C . calculate the total gradient, instead of just the C . intermolecular one C-RDC write(*,*) "iflaggr.eq.1 => performing tot grad" c inter = .false. c call gcart(grad) c if (pme) then c call pme_direct(edir1,eself1,ecoul1,grad) c call pme_recip(erec1,grad) c call pme_direct(edir1,eself1,ecoul1,gradpme) c call pme_recip(erec1,gradpme) c call pme_erec(erec,gradpme) c do i=1,natoms c i3 = 3*i c i2 = i3-1 c i1 = i2-1 c grad(i1) = grad(i1) c grad(i2) = grad(i2) c grad(i3) = grad(i3) c enddo C . add PME energies to total energy c elr = edir + eself + erec - ecoul c etot = etot + elr C . add PME contribution to heat of formation c eheat = eheat + 23.061*elr c endif c call setopt c call dograd(grad,gnorm,ierror) c if (ierror.ne.0) goto 1000 c endif C-RDC write(*,*) "edir1 eself1 ecoul1 erec1 ", edir1, eself1, ecoul1, erec1 do i=1,3*natoms !#ifdef MPI ! gradroar(i) = grade(i) / nproc !#else gradroar(i) = grade(i) !#endif enddo cccccccccccccccccccccccccccccccccccccccccccccccccc c call dovir(vir, gradroar) cccccccccccccccccccccccccccccccccccccccccccccccccc IF(prtvec) then IF (STAND)THEN call setprtvec NORBS = IORBPT(2) - 1 CALL WRTVEC(NORBS,EVEC,EVAL) elseif (frozenmc) then call setprtvec call wrtvecsubf endif ENDIF if (dos) call mkdos if(.not.CART)then call wrtint if(index(keywrd,'PRTCOORDS').ne.0) call wrtxyz else call wrtxyz endif if(GRADIENT) call wrtgrd(grade) call wrtchg c IF (pole) call dipoleh2o ! GM: 2010-10-31 begin ! add: call for dipole moment calculation if (pole) then call dipole endif ! GM: 2010-10-31 end IF (zgen) THEN call zmake call wrtint ENDIF IF (rotat.and..not.rotat2) call rotate 1000 return end