C C $Id: energy.F,v 1.8 1998/07/16 16:39:43 jjv5 Exp $ C C------------------------------------------------------------------------ SUBROUTINE ENERGY(newsub,EHEAT1,IERROR) C C DRIVER ROUTINE TO DETERMINE HEAT OF FORMATION EHEAT1. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" #include "constants.h" C #ifdef MPI #include "mpif.h" integer ijk #endif LOGICAL FIRST,newsub DATA FIRST /.TRUE./ SAVE FIRST,HFATMS real*8 dedir(3),derec(3),deelecsr(3) C C DENSITY MATRIX WILL GET INITIALIZED TO DIAGONAL FORM ON FIRST CALL. C IF(FIRST)THEN INIT = 0 ELSE INIT = 1 ENDIF C C C . UPDATE SUBSYSTEMS IF AUTOMATIC SUBSETTING IS BEING USED. C if (newsub) then C . avdv: gensub should be called before the call to energy C . because of the nsubmove option for the MC runs c CALL GENSUB(IERROR) c IF(IERROR.NE.0) RETURN c C . if a new subsetting is done, the frozen density matrix C . approximation doesn't work C . C . overrule when a frozen calculation is performed using dmx files if (.not.fdmx) then do i=1,nsub notfroz(i) = .true. enddo endif C C . PAIRLIST SHOULD BE UPDATED FOR D&C CALCULATION WITH AUTOMATIC C . SUBSETTING, OR FOR ANY CALCULATION IN WHICH THE CUTBOND KEYWORD C . HAS BEEN SPECIFIED, OR FOR THE VERY FIRST CALL TO ENERGY. C IF(autosub.OR. . INDEX(KEYWRD,'CUTBOND').NE.0.OR.FIRST)THEN C C . DO PRE-PAIRLIST SETUP. C CALL SPROC1(IERROR) IF(IERROR.NE.0) RETURN C C-RDC WRITE(IOUT,'(/" UPDATING PAIRLIST")') C-RDC c IF(SCREEN) WRITE(ISCR,'(/" UPDATING PAIRLIST")') C C . GENERATE PAIRLIST. C CALL BPAIR(IERROR) IF(IERROR.NE.0)THEN WRITE(IOUT,'(/" MAXIMUM ALLOWED NUMBER OF BONDED ATOM", . " PAIRS EXCEEDED AT ATOM ",I5, . /" INCREASE MBPAIR PARAMETER IN divcon.dim", . " AND RECOMPILE")') IERROR RETURN ENDIF C C ASSIGN POINTERS TO GLOBAL MATRICES (SOME DEPEND ON PAIRLIST). C CALL GLBPNT(IERROR) IF(IERROR.NE.0) RETURN C . THE DENSITY MATRIX SHOULD BE RE-INITIALIZED WHENEVER THE PAIRLIST C . IS UPDATED BECAUSE ITS STRUCTURE AND SIZE DEPENDS UPON THE C . PAIRLIST. C CALL INITP(INIT,ierror) if (ierror.ne.0) return C C C . CREATE ARRAYS THAT ALLOW QUICK LOCATION OF PAIRLIST ADDESSES C . OF SPECIFIED ATOM PAIRS. C CALL IJMAKE C C . GET NUMBERS OF SUBSYSTEMS SHARED BY EACH ATOM PAIR. C CALL SPROC2 endif C . endif (autosub, cutbond) endif C C CALCULATE AND STORE INTERATOMIC DISTANCES IF THIS IS NOT A C DIRECT CALCULATION. C IF(mcsim.or.(.not.DIRECT))THEN CALL RCALC(IERROR) IF(IERROR.NE.0) return ENDIF C C ON FIRST CALL, COMPUTE HEAT OF FORMATION OF SEPARATED ATOMS. ALSO C SET GRADIENT NORM TO ZERO SO THAT DEFAULT CONVERGENCE CRITERIA C WILL BE USED IN DOSCF. RETURN WITH IERROR=-1 IF USER HAS SPECIFIED C THE 'TESTRUN' KEYWORD. C IF(FIRST)THEN IF(INDEX(KEYWRD,'TESTRUN').NE.0)THEN C-RDC WRITE(IOUT,'(/" STOPPING TESTRUN CALCULATION")') IERROR = -1 RETURN ENDIF FIRST = .FALSE. HFATMS = 0.0D0 DO 100 I=1,NATOMS IAI = IATNUM(I) IF(IAI.EQ.0) GO TO 100 HFATMS = HFATMS + HFATM(IAI) - eV2kcal*EEATM(IAI) 100 CONTINUE GNORM = 0.0D0 ENDIF C C GET SCF ELECTRONIC ENERGY. C if (pmeqm) then c COMPUTE EWALD PAIR MATRIX do ij=1,natoms*natoms PHIPME(ij) = 0.0d0 DXPHIPME(ij) = 0.0d0 DYPHIPME(ij) = 0.0d0 DZPHIPME(ij) = 0.0d0 PHISR(ij) = 0.0d0 DXPHISR(ij) = 0.0d0 DYPHISR(ij) = 0.0d0 DZPHISR(ij) = 0.0d0 enddo call pme_setup c ij=0 #ifdef MPI DO I=myid+1,NATOMS,nproc #else do i=1,natoms #endif call pme_qm_direct(edir,eself,eclmb,dedir,i,i) call pme_qm_recip(erec,derec,i,i) ij=(i-1)*natoms+i PHIPME(ij) = erec + edir DXPHIPME(ij) = derec(1)+dedir(1) DYPHIPME(ij) = derec(2)+dedir(2) DZPHIPME(ij) = derec(3)+dedir(3) do j=1,i-1 call pme_qm_direct(edir,eself,eclmb,dedir,i,j) call pme_qm_recip(erec,derec,i,j) call elecsr(eelecsr,deelecsr,i,j) ij=(i-1)*natoms+j ji=(j-1)*natoms +i PHIPME(ij) = erec + edir PHISR(ij)=eelecsr PHIPME(ji)= PHIPME(ij) PHISR(ji)=PHISR(ij) DXPHIPME(ij) = derec(1)+dedir(1) DYPHIPME(ij) = derec(2)+dedir(2) DZPHIPME(ij) = derec(3)+dedir(3) DXPHIPME(ji) = -DXPHIPME(ij) DYPHIPME(ji) = -DYPHIPME(ij) DZPHIPME(ji) = -DZPHIPME(ij) DXPHISR(ij) = deelecsr(1) DYPHISR(ij) = deelecsr(2) DZPHISR(ij) = deelecsr(3) DXPHISR(ji) = -DXPHISR(ij) DYPHISR(ji) = -DYPHISR(ij) DZPHISR(ji) = -DZPHISR(ij) c write(59,*) edir , erec,PHIPME(ij) enddo enddo #ifdef MPI c IIMAX = IIMAT(NATOMS+1)-1 c IPMAX = IP1(NATOMS+1) c IJMAX = IJMAT(IPMAX)-1 call mpi_allreduce(PHIPME,TMP_MPI,natoms*natoms, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,natoms*natoms phipme(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(DXPHIPME,TMP_MPI,natoms*natoms, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,natoms*natoms dxphipme(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(DYPHIPME,TMP_MPI,natoms*natoms, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,natoms*natoms dyphipme(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(DZPHIPME,TMP_MPI,natoms*natoms, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,natoms*natoms dzphipme(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(PHISR,TMP_MPI,natoms*natoms, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,natoms*natoms phisr(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(DXPHISR,TMP_MPI,natoms*natoms, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,natoms*natoms dxphisr(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(DYPHISR,TMP_MPI,natoms*natoms, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,natoms*natoms dyphisr(ijk) = TMP_MPI(ijk) enddo call mpi_allreduce(DZPHISR,TMP_MPI,natoms*natoms, + MPI_DOUBLE_PRECISION,MPI_SUM,commsebomd,ier) do ijk = 1,natoms*natoms dzphisr(ijk) = TMP_MPI(ijk) enddo #endif c do i=1,natoms c do j=1,i-1 c call elecsr(eelecsr,deelecsr,i,j) c ij=(i-1)*natoms+j c ji=(j-1)*natoms +i c PHISR(ij)=eelecsr c PHISR(ji)=PHISR(ij) c DXPHISR(ij) = deelecsr(1) c DYPHISR(ij) = deelecsr(2) c DZPHISR(ij) = deelecsr(3) c DXPHISR(ji) = -DXPHISR(ij) c DYPHISR(ji) = -DYPHISR(ij) c DZPHISR(ji) = -DZPHISR(ij) c enddo c enddo endif IIMAX = IIMAT(NATOMS+1)-1 IPMAX = IP1(NATOMS+1) IJMAX = IJMAT(IPMAX)-1 nnorbs=0 do i=1,natoms nNORBS = nNORBS + NATORB(IATNUM(I)) enddo CALL DOSCF(EELECT,ECORE,ETOT,IERROR) C reject this configuration in case of MC-simulation and ierror.ne.0 if (ierror.ne.0) return C C COMPUTE MOLECULAR MECHANICS CORRECTION TO PEPTIDE TORSIONAL BARRIER C IF WE'VE STORED PEPTIDE BONDS. WILL BE IN eV. C 150 EPEP = 0.0D0 ! IF(N2PEP.GT.0)THEN ! DO 200 IPEP=1,N2PEP ! CALL DIHEDR(XYZ,IABC(1,IPEP),IABC(2,IPEP), ! . IABC(3,IPEP),IABC(4,IPEP),DIH) ! EIABC(IPEP) = EPEPMX*SIN(DIH)**2 ! EPEP = EPEP + EIABC(IPEP) ! 200 CONTINUE ! ENDIF C C HEAT OF FORMATION IN KCAL/MOL: C EHEAT1 = eV2kcal*(ETOT+EPEP) + HFATMS NSCF = NSCF + 1 RETURN END