hdiff output

r31689/md.f90 2017-01-21 10:35:55.302250000 +0000 r31688/md.f90 2017-01-21 10:35:55.582250000 +0000
 23:   !  23:   ! 
 24:   IMPLICIT NONE 24:   IMPLICIT NONE
 25:   ! 25:   !
 26:   ! Parameters set in the 'keywords' routine are to be saved! 26:   ! Parameters set in the 'keywords' routine are to be saved!
 27:   LOGICAL, SAVE :: MDT 27:   LOGICAL, SAVE :: MDT
 28:   INTEGER, SAVE :: MD_NWAIT, MD_NFREQ, MD_UNIT, MD_NSTEPS 28:   INTEGER, SAVE :: MD_NWAIT, MD_NFREQ, MD_UNIT, MD_NSTEPS
 29:   DOUBLE PRECISION, SAVE :: MD_GAMMA, MD_TSTEP 29:   DOUBLE PRECISION, SAVE :: MD_GAMMA, MD_TSTEP
 30:   ! 30:   !
 31:   ! MD-specific variables 31:   ! MD-specific variables
 32:   DOUBLE PRECISION :: EPOT, EKIN 32:   DOUBLE PRECISION :: EPOT, EKIN
 33:   DOUBLE PRECISION, ALLOCATABLE :: ACC(:), VEL(:)  33:   DOUBLE PRECISION, ALLOCATABLE :: FORCE(:), VEL(:) 
 34:   ! 34:   !
 35:   DOUBLE PRECISION, PARAMETER :: PI=3.141592652D0 35:   DOUBLE PRECISION, PARAMETER :: PI=3.141592652D0
 36:   ! 36:   !
 37: CONTAINS 37: CONTAINS
 38:   ! 38:   !
 39:   !============================================================ 39:   !============================================================
 40:   ! 40:   !
 41:   SUBROUTINE MDRUN(POS,NSTEPS,TEMP,DUMP) 41:   SUBROUTINE MDRUN(POS,NSTEPS,TEMP,DUMP)
 42:     ! 42:     !
 43:     IMPLICIT NONE 43:     IMPLICIT NONE
 46:     INTEGER, INTENT(IN) :: NSTEPS 46:     INTEGER, INTENT(IN) :: NSTEPS
 47:     DOUBLE PRECISION, INTENT(IN) :: TEMP ! k_B*T in units of energy 47:     DOUBLE PRECISION, INTENT(IN) :: TEMP ! k_B*T in units of energy
 48:     DOUBLE PRECISION, INTENT(INOUT) :: POS(3*NATOMS) 48:     DOUBLE PRECISION, INTENT(INOUT) :: POS(3*NATOMS)
 49:     ! 49:     !
 50:     INTEGER :: I, J, K 50:     INTEGER :: I, J, K
 51:     DOUBLE PRECISION :: TIME 51:     DOUBLE PRECISION :: TIME
 52:     ! 52:     !
 53:     IF (DUMP) WRITE(MYUNIT, '(A)') "mdrun> Start of MD run." 53:     IF (DUMP) WRITE(MYUNIT, '(A)') "mdrun> Start of MD run."
 54:         54:        
 55:     ! 55:     !
 56:     ALLOCATE(ACC(3*NATOMS), VEL(3*NATOMS)) 56:     ALLOCATE(FORCE(3*NATOMS), VEL(3*NATOMS))
 57:     ! 57:     !
 58:     CALL INIVEL(POS, TEMP)  58:     CALL INIVEL(POS, TEMP) 
 59:     CALL POTENTIAL(POS,ACC,EPOT,.TRUE.,.FALSE.) ! get gradient 59:     CALL POTENTIAL(POS,FORCE,EPOT,.TRUE.,.FALSE.) ! get gradient
 60:     K=0 60:     FORCE(:) = -FORCE(:) ! get force    
 61:     DO I=1,NATOMS 
 62:        DO J=1,3 
 63:           K=K+1 
 64:           ACC(K) = -ACC(K)/ATMASS(I) ! get ACC 
 65:        ENDDO 
 66:     ENDDO 
 67:     ! 61:     !
 68:     I=0 62:     I=0
 69:     IF (DUMP) WRITE(MYUNIT, '(A,I8,A,F18.8,A,F14.8)') & 63:     IF (DUMP) WRITE(MYUNIT, '(A,I8,A,F18.8,A,F14.8)') &
 70:          'mdrun> Step ', I,' EPOT= ',EPOT,' EKIN= ',EKIN 64:          'mdrun> Step ', I,' EPOT= ',EPOT,' EKIN= ',EKIN
 71:     ! 65:     !
 72:     DO I=1,NSTEPS 66:     DO I=1,NSTEPS
 73:        ! 67:        !
 74:        CALL MDSTEP(POS, TEMP) 68:        CALL MDSTEP(POS, TEMP)
 75:        ! 69:        !
 76:        IF(DUMP) THEN ! Dump data to files 70:        IF(DUMP) THEN ! Dump data to files
 92:              DO J=1,NATOMS 86:              DO J=1,NATOMS
 93:                 WRITE(DUMPXYZUNIT(1),'(A,3(1X,F10.5))') & 87:                 WRITE(DUMPXYZUNIT(1),'(A,3(1X,F10.5))') &
 94:                      'X ', (POS(3*(J-1)+K), K=1,3) 88:                      'X ', (POS(3*(J-1)+K), K=1,3)
 95:              ENDDO              89:              ENDDO             
 96:           ENDIF 90:           ENDIF
 97:           ! 91:           !
 98:        ENDIF 92:        ENDIF
 99:        ! 93:        !
100:     ENDDO 94:     ENDDO
101:     ! 95:     !
102:     DEALLOCATE(ACC, VEL) 96:     DEALLOCATE(FORCE, VEL)
103:     ! 97:     !
104:     CALL MYCPU_TIME(TIME) 98:     CALL MYCPU_TIME(TIME)
105:     IF (DUMP) WRITE(MYUNIT, '(A,F11.1)') & 99:     IF (DUMP) WRITE(MYUNIT, '(A,F11.1)') &
106:          "mdrun> Finished MD run t= ", TIME-TSTART100:          "mdrun> Finished MD run t= ", TIME-TSTART
107:     !101:     !
108:     RETURN102:     RETURN
109:     !103:     !
110:   END SUBROUTINE MDRUN104:   END SUBROUTINE MDRUN
111:   !105:   !
112:   !============================================================106:   !============================================================
123:          DPRAND,R1,R2,NR,NR2(3*NATOMS)117:          DPRAND,R1,R2,NR,NR2(3*NATOMS)
124:     !118:     !
125:     HDT=0.5D0*MD_TSTEP119:     HDT=0.5D0*MD_TSTEP
126:     GFRIC=1.0D0-MD_GAMMA*HDT120:     GFRIC=1.0D0-MD_GAMMA*HDT
127:     NOISE(:) = DSQRT(TEMP*MD_GAMMA*MD_TSTEP/ATMASS(:))121:     NOISE(:) = DSQRT(TEMP*MD_GAMMA*MD_TSTEP/ATMASS(:))
128:     !122:     !
129:     K=0123:     K=0
130:     DO I=1,NATOMS124:     DO I=1,NATOMS
131:        DO J=1,3125:        DO J=1,3
132:           K=K+1126:           K=K+1
133:           ! Generate normal random numbers127:           ! Generate norman random numbers
134:           R1=DSQRT(-2.0D0*LOG(DPRAND()))128:           R1=DSQRT(-2.0D0*LOG(DPRAND()))
135:           R2=2.0*PI*DPRAND()129:           R2=2.0*PI*DPRAND()
136:           NR=R1*COS(R2)130:           NR=R1*COS(R2)
137:           NR2(K)=R1*SIN(R2)          131:           NR2(K)=R1*SIN(R2)          
138:           ! Update mid-step velocity132:           ! Update mid-step velocity
139:           VEL(K) = GFRIC*VEL(K) + ACC(K)*HDT + NR*NOISE(I)133:           VEL(K) = GFRIC*VEL(K) + FORCE(K)*HDT/ATMASS(I) &
 134:                + NR*NOISE(I)
140:           ! Update full-step position135:           ! Update full-step position
141:           POS(K) = POS(K) + VEL(K)*MD_TSTEP136:           POS(K) = POS(K) + VEL(K)*MD_TSTEP
142:        ENDDO137:        ENDDO
143:     ENDDO138:     ENDDO
144:     !139:     !
145:     CALL POTENTIAL(POS,ACC,EPOT,.TRUE.,.FALSE.) ! GET GRADIENT140:     CALL POTENTIAL(POS,FORCE,EPOT,.TRUE.,.FALSE.) ! GET GRADIENT
 141:     FORCE(:)= -FORCE(:)
146:     !142:     !
147:     K=0143:     K=0
148:     EKIN=0.0D0144:     EKIN=0.0D0
149:     DO I=1,NATOMS145:     DO I=1,NATOMS
150:        DO J=1,3146:        DO J=1,3
151:           K=K+1147:           K=K+1
152:           ACC(K)= -ACC(K)/ATMASS(I) ! get ACC 
153:           ! Update full-step velocity148:           ! Update full-step velocity
154:           VEL(K) = GFRIC*VEL(K) + ACC(K)*HDT + NR2(K)*NOISE(I)149:           VEL(K) = GFRIC*VEL(K) + FORCE(K)*HDT/ATMASS(I) &
 150:                + NR2(K)*NOISE(I)
155:           EKIN = EKIN + ATMASS(I)*VEL(K)*VEL(K)151:           EKIN = EKIN + ATMASS(I)*VEL(K)*VEL(K)
156:        ENDDO152:        ENDDO
157:     ENDDO153:     ENDDO
158:     EKIN=0.5D0*EKIN154:     EKIN=0.5D0*EKIN
159:     !155:     !
160:     RETURN156:     RETURN
161:     !157:     !
162:   END SUBROUTINE MDSTEP158:   END SUBROUTINE MDSTEP
163:   !159:   !
164:   !============================================================160:   !============================================================


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0