hdiff output

r31117/keywords.f 2016-09-05 18:30:14.249076162 +0100 r31116/keywords.f 2016-09-05 18:30:15.005086504 +0100
 42:       USE MBPOLMOD, ONLY: MBPOLINIT 42:       USE MBPOLMOD, ONLY: MBPOLINIT
 43:       USE SWMOD, ONLY: SWINIT, MWINIT 43:       USE SWMOD, ONLY: SWINIT, MWINIT
 44:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS, 44:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS,
 45:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA 45:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA
 46:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 46:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 47:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 47:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 48:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,  48:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 
 49:      &     PARSE_MLJ_PARAMS 49:      &     PARSE_MLJ_PARAMS
 50:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT 50:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT
 51:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE 51:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE
 52:       USE MD, ONLY : MDT, MD_TSTEP, MD_GAMMA, MD_NWAIT, MD_NFREQ, MD_NSTEPS 
 53:        
 54:       IMPLICIT NONE 52:       IMPLICIT NONE
 55:  53: 
 56:       DOUBLE PRECISION, ALLOCATABLE :: MLPMEAN(:) 54:       DOUBLE PRECISION, ALLOCATABLE :: MLPMEAN(:)
 57:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, LAST, IX, J1, JP, NPCOUNT, NDUMMY, INDEX, J2, J3, J4 55:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, LAST, IX, J1, JP, NPCOUNT, NDUMMY, INDEX, J2, J3, J4
 58:       INTEGER DATA_UNIT 56:       INTEGER DATA_UNIT
 59:       INTEGER MOVABLEATOMINDEX 57:       INTEGER MOVABLEATOMINDEX
 60:       LOGICAL CAT, YESNO, PERMFILE, CONFILE 58:       LOGICAL CAT, YESNO, PERMFILE, CONFILE
 61:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR, 59:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR,
 62:      &                NERROR, ECHO, LAST, CAT 60:      &                NERROR, ECHO, LAST, CAT
 63:        DOUBLE PRECISION XX, ROH, ROM, WTHETA  61:        DOUBLE PRECISION XX, ROH, ROM, WTHETA 
550:       SAWTOOTH_TFAC=1.0D0548:       SAWTOOTH_TFAC=1.0D0
551:       SAWTOOTH_SFAC=1.0D0549:       SAWTOOTH_SFAC=1.0D0
552:       SAWTOOTH_SFAC2=1.0D0550:       SAWTOOTH_SFAC2=1.0D0
553: 551: 
554: !     ds656> Generalised LJ with Yukawa552: !     ds656> Generalised LJ with Yukawa
555:       GLJY = .FALSE.553:       GLJY = .FALSE.
556:       GLJ_EXP = 6554:       GLJ_EXP = 6
557:       YUK_A = 0.0D0555:       YUK_A = 0.0D0
558:       YUK_XI = 1.0D0556:       YUK_XI = 1.0D0
559:       557:       
560: !     ds656> Molecular Dynamics 
561:       MDT = .FALSE. 
562:       MD_NSTEPS = 0 
563:       MD_NWAIT = 1 
564:       MD_NFREQ = 100 
565:       MD_GAMMA = 0.1 
566:       MD_TSTEP = 0.01 
567:  
568:       CHRMMT=.FALSE.558:       CHRMMT=.FALSE.
569:       CHARMMTYPE=1559:       CHARMMTYPE=1
570:       CHARMMDFTBT=.FALSE.560:       CHARMMDFTBT=.FALSE.
571:       ACESOLV=.FALSE.561:       ACESOLV=.FALSE.
572:       ACEUPSTEP=50562:       ACEUPSTEP=50
573:       CHRIGIDTRANST=.FALSE.563:       CHRIGIDTRANST=.FALSE.
574:       CHRIGIDROTT=.FALSE.564:       CHRIGIDROTT=.FALSE.
575:       CHNMIN=0.D0565:       CHNMIN=0.D0
576:       CHNMAX=HUGE(1.0D0)566:       CHNMAX=HUGE(1.0D0)
577:       CHMDT=.FALSE.567:       CHMDT=.FALSE.
2119:         ELSE2109:         ELSE
2120:             CALL MMEINITWRAPPER(TRIM(ADJUSTL(PRMTOP))//C_NULL_CHAR, IGB, SALTCON, RGBMAX, CUT)2110:             CALL MMEINITWRAPPER(TRIM(ADJUSTL(PRMTOP))//C_NULL_CHAR, IGB, SALTCON, RGBMAX, CUT)
2121:         END IF2111:         END IF
2122: 2112: 
2123: !2113: !
2124: ! The maximum number of constraints to use in the constrained potential.2114: ! The maximum number of constraints to use in the constrained potential.
2125: ! The deafult is 4.2115: ! The deafult is 4.
2126: !2116: !
2127:       ELSE IF (WORD.EQ.'MAXCON') THEN2117:       ELSE IF (WORD.EQ.'MAXCON') THEN
2128:          CALL READI(MAXCONUSE)2118:          CALL READI(MAXCONUSE)
2129:  
2130: ! ds656> Molecular Dynamics 
2131:       ELSE IF (WORD.EQ.'MD') THEN 
2132:          MDT=.TRUE. 
2133:          CALL READI(MD_NSTEPS) 
2134:          IF(NITEMS .GT. 2) THEN 
2135:             CALL READI(MD_NWAIT) 
2136:             CALL READI(MD_NFREQ) 
2137:          ELSE 
2138:             MD_NWAIT=MD_NSTEPS 
2139:             MD_NFREQ=1 
2140:          ENDIF 
2141:           
2142:       ELSE IF (WORD.EQ.'MDPARAMS') THEN 
2143:          CALL READF(MD_TSTEP) 
2144:          CALL READF(MD_GAMMA) 
2145: !2119: !
2146: ! MK ion trap2120: ! MK ion trap
2147: !2121: !
2148:       ELSE IF (WORD.EQ.'MKTRAP') THEN2122:       ELSE IF (WORD.EQ.'MKTRAP') THEN
2149:          MKTRAPT=.TRUE.2123:          MKTRAPT=.TRUE.
2150: !2124: !
2151: ! Three layer neural network (multilayer perceptron) with2125: ! Three layer neural network (multilayer perceptron) with
2152: ! MLPIN inputs (columns per data item)2126: ! MLPIN inputs (columns per data item)
2153: ! MLPOUT outputs2127: ! MLPOUT outputs
2154: ! MLPHIDDEN hidden nodes2128: ! MLPHIDDEN hidden nodes


r31117/main.F 2016-09-05 18:30:14.501079608 +0100 r31116/main.F 2016-09-05 18:30:15.257089946 +0100
 27:       USE MODAMBER 27:       USE MODAMBER
 28:       USE MODAMBER9, only : AMBFINALIO_NODE,MDCRD_UNIT,MDINFO_UNIT,AMBPDB_UNIT, ATMASS1 28:       USE MODAMBER9, only : AMBFINALIO_NODE,MDCRD_UNIT,MDINFO_UNIT,AMBPDB_UNIT, ATMASS1
 29:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_MASSES 29:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_MASSES
 30:       USE MODCHARMM 30:       USE MODCHARMM
 31:       USE PORFUNCS 31:       USE PORFUNCS
 32:       USE TWIST_MOD 32:       USE TWIST_MOD
 33:       USE HOMOREFMOD 33:       USE HOMOREFMOD
 34:       USE GENRIGID, only : RIGIDINIT, GENRIGID_READ_FROM_FILE, DEGFREEDOMS 34:       USE GENRIGID, only : RIGIDINIT, GENRIGID_READ_FROM_FILE, DEGFREEDOMS
 35:       USE MULTIPOT, only: MULTIPOT_INITIALISE 35:       USE MULTIPOT, only: MULTIPOT_INITIALISE
 36:       USE GAUSS_MOD, ONLY: KEGEN 36:       USE GAUSS_MOD, ONLY: KEGEN
 37:       USE MD, ONLY : MDT, MD_NSTEPS, MD_NWAIT, MDRUN 
 38:  
 39:       IMPLICIT NONE 37:       IMPLICIT NONE
 40:       !EXTERNAL READ_CMD_ARGS 38:       !EXTERNAL READ_CMD_ARGS
 41: #ifdef MPI 39: #ifdef MPI
 42:       INCLUDE 'mpif.h' 40:       INCLUDE 'mpif.h'
 43: #endif 41: #endif
 44:       INTEGER J1,J2, JP, MPIERR, NDUMMY3,NPTOTAL,VERSIONTEMP,GETUNIT,LUNIT 42:       INTEGER J1,J2, JP, MPIERR, NDUMMY3,NPTOTAL,VERSIONTEMP,GETUNIT,LUNIT
 45:       DOUBLE PRECISION, ALLOCATABLE :: SCREENC(:) 43:       DOUBLE PRECISION, ALLOCATABLE :: SCREENC(:)
 46:       DOUBLE PRECISION POTEL, DUMMY, INTFREEZETOLSAVE, RMAT(3,3), DUMMY2, DIST2, LINTCONSTRAINTTOL 44:       DOUBLE PRECISION POTEL, DUMMY, INTFREEZETOLSAVE, RMAT(3,3), DUMMY2, DIST2, LINTCONSTRAINTTOL
 47:       DOUBLE PRECISION, ALLOCATABLE :: TMPCOORDS(:), ENDCOORDS(:,:) 45:       DOUBLE PRECISION, ALLOCATABLE :: TMPCOORDS(:), ENDCOORDS(:,:)
 48:       INTEGER, ALLOCATABLE :: NDUMMY(:), NDUMMY2(:,:) 46:       INTEGER, ALLOCATABLE :: NDUMMY(:), NDUMMY2(:,:)
373:       FF(1:NSAVE)=0 ! to prevent reading from uninitialised memorY371:       FF(1:NSAVE)=0 ! to prevent reading from uninitialised memorY
374:       VATO(1:NATOMSALLOC,1:NPAR)=0.0D0 ! to prevent reading from uninitialised memory372:       VATO(1:NATOMSALLOC,1:NPAR)=0.0D0 ! to prevent reading from uninitialised memory
375:       ALLOCATE(ESAVE(NTAB,NPAR),XINSAVE(NTAB,NPAR))373:       ALLOCATE(ESAVE(NTAB,NPAR),XINSAVE(NTAB,NPAR))
376:       ALLOCATE(VEC(NVEC))374:       ALLOCATE(VEC(NVEC))
377: 375: 
378: !     IF (SYMMETRIZE.AND.(.NOT.CENT)) THEN376: !     IF (SYMMETRIZE.AND.(.NOT.CENT)) THEN
379: !        PRINT '(A)','Probable input error - SYMMETRIZE true but CENT false'377: !        PRINT '(A)','Probable input error - SYMMETRIZE true but CENT false'
380: !        STOP378: !        STOP
381: !     ENDIF379: !     ENDIF
382: 380: 
383:       IF (DUMPT .OR. DUMP_MARKOV .OR. (MD_NWAIT.LE.MD_NSTEPS)) THEN381:       IF (DUMPT .OR. DUMP_MARKOV) THEN
384:          IF (CHRMMT) THEN382:          IF (CHRMMT) THEN
385:             OPEN(UNIT=719,FILE='dump.crd',STATUS='UNKNOWN')383:             OPEN(UNIT=719,FILE='dump.crd',STATUS='UNKNOWN')
386:             OPEN(UNIT=720,FILE='dump.pdb',STATUS='UNKNOWN')384:             OPEN(UNIT=720,FILE='dump.pdb',STATUS='UNKNOWN')
387:          ENDIF385:          ENDIF
388: !386: !
389: ! dump.1.xyz is partly filled with control characters for no apparent reason.387: ! dump.1.xyz is partly filled with control characters for no apparent reason.
390: ! Suspect a compiler or mpi bug? DJW388: ! Suspect a compiler or mpi bug? DJW
391: ! dump.1.xyz is fine with debug compilation, so it is a compiler bug!389: ! dump.1.xyz is fine with debug compilation, so it is a compiler bug!
392: !390: !
393:          ALLOCATE(DUMPXYZUNIT(NPAR),DUMPVUNIT(NPAR))391:          ALLOCATE(DUMPXYZUNIT(NPAR),DUMPVUNIT(NPAR))
581:          NPCALL_QMIN(J1)=0579:          NPCALL_QMIN(J1)=0
582:       ENDDO580:       ENDDO
583: 581: 
584:       CALL INITIALIZATIONS()582:       CALL INITIALIZATIONS()
585: 583: 
586: ! vr274>  call test routines.584: ! vr274>  call test routines.
587: !         By default this routine does nothing. In the testsuite, it is replaced by585: !         By default this routine does nothing. In the testsuite, it is replaced by
588: !         a customized function which can run checks whether everything runs fine.586: !         a customized function which can run checks whether everything runs fine.
589:       CALL RUN_TESTS_AFTER_INIT()587:       CALL RUN_TESTS_AFTER_INIT()
590: 588: 
591:       IF (GENALT) THEN ! mo361> Run GA589: !mo361> Run GA
592:          !NRUNS=0590:       IF (GENALT) THEN
 591:          NRUNS=0
593:          CALL MYGA_RUN()592:          CALL MYGA_RUN()
594:       ELSEIF (MULTIPERMT) THEN  ! ds656> Span label permutations593:       ENDIF
 594: !ds656> Span label permutations
 595:       IF (MULTIPERMT) THEN
 596:          NRUNS=0
595:          CALL MULTIPERM()597:          CALL MULTIPERM()
596:       ELSEIF (MDT) THEN          ! ds656> Molecular Dynamics 
597:          CALL MDRUN(COORDS,MD_NSTEPS,TEMP(1),.TRUE.) 
598:       ELSEIF ((NRUNS.GT.0).OR.PTMC.OR.BSPT) THEN 
599:          CALL MCRUNS(SCREENC) 
600:       ENDIF598:       ENDIF
 599: 
 600:       IF ((NRUNS.GT.0).OR.PTMC.OR.BSPT) CALL MCRUNS(SCREENC)
601: C     CALL SYSTEM('rm ssdump ssave >& /dev/null')601: C     CALL SYSTEM('rm ssdump ssave >& /dev/null')
602: 602: 
603:       IF (ALLOCATED(FIN)) DEALLOCATE(FIN)603:       IF (ALLOCATED(FIN)) DEALLOCATE(FIN)
604:       IF (ALLOCATED(XICOM)) DEALLOCATE(XICOM)604:       IF (ALLOCATED(XICOM)) DEALLOCATE(XICOM)
605:       IF (ALLOCATED(PCOM)) DEALLOCATE(PCOM)605:       IF (ALLOCATED(PCOM)) DEALLOCATE(PCOM)
606:       IF (ALLOCATED(GAUSSKK)) DEALLOCATE(GAUSSKK,GAUSSEE)606:       IF (ALLOCATED(GAUSSKK)) DEALLOCATE(GAUSSKK,GAUSSEE)
607: C     deallocate dynamic memory for AMBER607: C     deallocate dynamic memory for AMBER
608:       IF (AMBERT) CALL AMBER_DEALLOCATE_STACKS608:       IF (AMBERT) CALL AMBER_DEALLOCATE_STACKS
609:       IF (ALLOCATED(ANV)) DEALLOCATE(ANV)      609:       IF (ALLOCATED(ANV)) DEALLOCATE(ANV)      
610:       CALL FLUSH(MYUNIT)610:       CALL FLUSH(MYUNIT)


r31117/md.f90 2016-09-05 18:30:14.745082946 +0100 r31116/md.f90 2016-09-05 18:30:15.509093397 +0100
  1: !  GMIN: A program for finding global minima  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/md.f90' in revision 31116
  2: !  Copyright (C) 1999-2006 David J. Wales 
  3: !  This file is part of GMIN. 
  4: ! 
  5: !  GMIN is free software; you can redistribute it and/or modify 
  6: !  it under the terms of the GNU General Public License as published by 
  7: !  the Free Software Foundation; either version 2 of the License, or 
  8: !  (at your option) any later version. 
  9: ! 
 10: !  GMIN is distributed in the hope that it will be useful, 
 11: !  but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: !  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: !  GNU General Public License for more details. 
 14: ! 
 15: !  You should have received a copy of the GNU General Public License 
 16: !  along with this program; if not, write to the Free Software 
 17: !  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA 
 18: ! 
 19: MODULE MD 
 20:   ! 
 21:   USE COMMONS, ONLY : NATOMS, MYUNIT, ATMASS, PRTFRQ, DUMPXYZUNIT 
 22:   !  
 23:   IMPLICIT NONE 
 24:   ! 
 25:   ! Parameters set by keywords 
 26:   LOGICAL :: MDT 
 27:   INTEGER :: MD_NWAIT, MD_NFREQ, MD_UNIT, MD_NSTEPS 
 28:   DOUBLE PRECISION :: MD_GAMMA, MD_TSTEP 
 29:   ! 
 30:   ! MD-specific variables 
 31:   DOUBLE PRECISION :: EPOT, EKIN 
 32:   DOUBLE PRECISION, ALLOCATABLE :: FORCE(:), VEL(:)  
 33:   ! 
 34: CONTAINS 
 35:   ! 
 36:   !============================================================ 
 37:   ! 
 38:   SUBROUTINE MDRUN(POS,NSTEPS,TEMP,DUMP) 
 39:     ! 
 40:     IMPLICIT NONE 
 41:     ! 
 42:     LOGICAL, INTENT(IN) :: DUMP 
 43:     INTEGER, INTENT(IN) :: NSTEPS 
 44:     DOUBLE PRECISION, INTENT(IN) :: TEMP ! k_B*T in units of energy 
 45:     DOUBLE PRECISION, INTENT(INOUT) :: POS(3*NATOMS) 
 46:     ! 
 47:     INTEGER :: I, J, K 
 48:     ! 
 49:     IF (DUMP) WRITE(MYUNIT, '(A)') "mdrun> Start of MD run." 
 50:         
 51:     ! 
 52:     ALLOCATE(FORCE(3*NATOMS), VEL(3*NATOMS)) 
 53:     ! 
 54:     CALL INIVEL(POS, TEMP)  
 55:     CALL POTENTIAL(POS,FORCE,EPOT,.TRUE.,.FALSE.) ! get gradient 
 56:     FORCE(:) = -FORCE(:) ! get force     
 57:     ! 
 58:     I=0 
 59:     IF (DUMP) WRITE(MYUNIT, '(A,I8,A,F18.8,A,F14.8)') & 
 60:          'mdrun> Step ', I,' EPOT= ',EPOT,' EKIN= ',EKIN 
 61:     ! 
 62:     DO I=1,NSTEPS 
 63:        ! 
 64:        CALL MDSTEP(POS, TEMP) 
 65:        ! 
 66:        IF(DUMP) THEN ! Dump data to files 
 67:           ! 
 68:           IF(MOD(I,PRTFRQ).EQ.0) THEN ! Basic data 
 69:              WRITE(MYUNIT, '(A,I8,A,F18.8,A,F14.8)') & 
 70:                   'mdrun> Step ', I,' EPOT= ',EPOT,' EKIN= ',EKIN 
 71:           ENDIF 
 72:           ! 
 73:           IF(I > MD_NWAIT .AND. & 
 74:                MOD(I-MD_NWAIT,MD_NFREQ)==0) THEN ! XYZ dump 
 75:              WRITE(DUMPXYZUNIT(1),'(I4)') NATOMS 
 76:              WRITE(DUMPXYZUNIT(1),'(A,I9,A,F20.10)')  & 
 77:                   'MD step ',I,'  EPOT= ', EPOT 
 78:              DO J=1,NATOMS 
 79:                 WRITE(DUMPXYZUNIT(1),'(A,3(1X,F10.5))') & 
 80:                      'X ', (POS(3*(J-1)+K), K=1,3) 
 81:              ENDDO              
 82:           ENDIF 
 83:           ! 
 84:        ENDIF 
 85:        ! 
 86:     ENDDO 
 87:     ! 
 88:     DEALLOCATE(FORCE, VEL) 
 89:     ! 
 90:     IF (DUMP) WRITE(MYUNIT, '(A)') "mdrun> End of MD run." 
 91:     ! 
 92:     RETURN 
 93:     ! 
 94:   END SUBROUTINE MDRUN 
 95:   ! 
 96:   !============================================================ 
 97:   ! 
 98:   SUBROUTINE MDSTEP(POS, TEMP) 
 99:     ! 
100:     IMPLICIT NONE 
101:     ! 
102:     DOUBLE PRECISION, INTENT(IN) :: TEMP 
103:     DOUBLE PRECISION, INTENT(INOUT) :: POS(3*NATOMS) 
104:      
105:     INTEGER :: I,J,K 
106:     DOUBLE PRECISION :: HDT, HDT2, GFRIC, NOISE(NATOMS), & 
107:          DPRAND 
108:     ! 
109:     HDT=0.5D0*MD_TSTEP 
110:     HDT2=HDT*MD_TSTEP 
111:     GFRIC=1.0D0-MD_GAMMA*HDT 
112:     NOISE(:) = DSQRT(6.0D0*TEMP*ATMASS(:)*MD_GAMMA/MD_TSTEP) 
113:     ! 
114:     K=0 
115:     DO I=1,NATOMS 
116:        DO J=1,3 
117:           K=K+1 
118:           POS(K) = POS(K) + VEL(K)*MD_TSTEP + FORCE(K)*HDT2/ATMASS(I) 
119:           VEL(K) = GFRIC*VEL(K) + FORCE(K)*HDT/ATMASS(I) 
120:        ENDDO 
121:     ENDDO 
122:     ! 
123:     CALL POTENTIAL(POS,FORCE,EPOT,.TRUE.,.FALSE.) ! GET GRADIENT 
124:     ! 
125:     K=0 
126:     EKIN=0.0D0 
127:     DO I=1,NATOMS 
128:        DO J=1,3 
129:           K=K+1 
130:           FORCE(K)= -FORCE(K) + (2.0D0*DPRAND()-1.0D0)*NOISE(I) 
131:           VEL(K) = GFRIC*VEL(K) + FORCE(K)*HDT/ATMASS(I) 
132:           EKIN = EKIN + ATMASS(I)*VEL(K)*VEL(K) 
133:        ENDDO 
134:     ENDDO 
135:     EKIN=0.5D0*EKIN 
136:     ! 
137:     RETURN 
138:     ! 
139:   END SUBROUTINE MDSTEP 
140:   ! 
141:   !============================================================ 
142:   ! 
143:   SUBROUTINE INIVEL(POS, TEMP) 
144:     ! 
145:     IMPLICIT NONE 
146:     ! 
147:     DOUBLE PRECISION, INTENT (IN) :: TEMP, POS(3*NATOMS) 
148:     ! 
149:     INTEGER :: I, J, K, L, IPIVOT(3) 
150:     DOUBLE PRECISION :: COM(3), COM_MOM(3), ANG_MOM(3), & 
151:          MOI(3,3), MOI_INV(3,3), DPRAND, SCALE, TMASS, & 
152:          P(3), R(3), W(3), R2, X(NATOMS,3) 
153:     DOUBLE PRECISION, PARAMETER :: PI=3.141592652D0 
154:     ! 
155:     COM(:) = 0.0D0 
156:     COM_MOM(:) = 0.0D0 
157:     TMASS = 0.0D0 
158:     ! 
159:     ! On first pass randomly generate velocities 
160:     K=0     
161:     DO I=1,NATOMS 
162:        TMASS = TMASS + ATMASS(I)  
163:        DO J=1,3 
164:           K=K+1 
165:           ! Pick random number from normal distribution with mean of 
166:           ! of zero and unit variance (See Allen and Tildesley p347) 
167:           VEL(K) = DSQRT(-2.0D0*LOG(DBLE(DPRAND())))* & 
168:                COS(2.0D0*PI*DBLE(DPRAND())) 
169:           ! 
170:           ! Will need to subtract net translation and rotation! 
171:           COM(J) = COM(J) + POS(K)*ATMASS(I) 
172:           COM_MOM(J) = COM_MOM(J) + VEL(K)*ATMASS(I) 
173:           ! 
174:        ENDDO 
175:     ENDDO 
176:     COM(:) = COM(:)/TMASS 
177:     ! 
178:     WRITE(MYUNIT,'(A,3(1x,F12.6))') & 
179:          'inivel> Centre of mass:', COM 
180:     WRITE(MYUNIT,'(A,3(1x,F12.6))') & 
181:          'inivel> Initial COM_MOM:', COM_MOM 
182:     ! On second pass substract net translation and determine rotation 
183:     ANG_MOM(:) = 0.0D0 ! L = Iw 
184:     MOI(:,:) = 0.0D0 
185:     K=0     
186:     DO I=1,NATOMS 
187:        R2=0.0D0 
188:        DO J=1,3 
189:           K=K+1 
190:           VEL(K) = VEL(K) - COM_MOM(J)/TMASS 
191:           R(J) = POS(K) - COM(J)  ! pos rel. to CoM 
192:           X(I,J) = R(J) 
193:           P(J) = VEL(K)*ATMASS(I) ! momentum component 
194:           R2 = R2 + R(J)*R(J) 
195:        ENDDO 
196:        ANG_MOM = ANG_MOM + CROSS(R,P) 
197:        DO J=1,3 
198:           MOI(J,J) = MOI(J,J) + (R2 - R(J)*R(J))*ATMASS(I) ! diagonal 
199:           DO L=J+1,3 ! off-diagonal 
200:              MOI(J,L) = MOI(J,L) - R(J)*R(L)*ATMASS(I) 
201:              MOI(L,J) = MOI(J,L) ! impose symmetry 
202:           ENDDO 
203:        ENDDO 
204:     ENDDO 
205:     ! 
206:     WRITE(MYUNIT,'(A,3(1X,F12.6))') & 
207:          'inivel> Initial ANG_MOM=', ANG_MOM 
208:     ! 
209:     ! Invert MOI using routines from LAPACK  
210:     CALL DGETRF(3,3,MOI,3,IPIVOT,I) ! MOI is modified!!!    
211:     CALL DGETRI(3,MOI,3,IPIVOT,R,3,I) 
212:     ! 
213:     IF(I /= 0) THEN 
214:        WRITE(MYUNIT,'(A,I1)') & 
215:             'md> Failed to invert inertia tensor! Error:', I 
216:        STOP 
217:     ENDIF 
218:     ! 
219:     ! Compute net angular velocity and store in W 
220:     DO I=1,3 
221:        W(I)=DOT_PRODUCT(MOI(I,1:3), ANG_MOM(1:3))  
222:     ENDDO 
223:     ! 
224:     ! On third pass eliminate rotation and compute kinetic energy 
225:     EKIN=0 
226:     K=0 
227:     ANG_MOM(:) = 0.0D0 ! L = Iw 
228:     COM_MOM(:) = 0.0D0 
229:     DO I=1,NATOMS 
230:        R=CROSS(W,X(I,1:3)) 
231:        DO J=1,3 
232:           K=K+1 
233:           VEL(K) = VEL(K) - R(J) 
234:           EKIN = EKIN + ATMASS(I)*VEL(K)*VEL(K) 
235:           P(J) = VEL(K)*ATMASS(I) ! momentum component 
236:           COM_MOM(J) = COM_MOM(J) + VEL(K)*ATMASS(I) 
237:        ENDDO 
238:        ANG_MOM = ANG_MOM + CROSS(X(I,1:3),P) 
239:        !write(myunit, *) ang_mom 
240:     ENDDO 
241:     ! 
242:     WRITE(MYUNIT,'(A,3(1x,F12.6))') & 
243:          'inivel> Final COM_MOM=', COM_MOM 
244:     WRITE(MYUNIT,'(A,3(1x,F12.6))') & 
245:          'inivel> Final ANG_MOM=', ANG_MOM 
246:     ! 
247:     SCALE=DSQRT(DBLE(K)*TEMP/EKIN) 
248:     ! 
249:     ! On fourth pass rescale velocities to correct temperature 
250:     K=0  
251:     EKIN=0.D0 
252:     DO I=1,NATOMS 
253:        DO J=1,3 
254:           K=K+1 
255:           VEL(K) = SCALE*VEL(K) 
256:           EKIN = EKIN + ATMASS(I)*VEL(K)*VEL(K) 
257:        ENDDO 
258:     ENDDO 
259:     EKIN=0.5D0*EKIN 
260:     ! 
261:     RETURN 
262:     ! 
263:   END SUBROUTINE INIVEL 
264:   ! 
265:   !============================================================ 
266:   ! 
267:   FUNCTION CROSS(A,B) 
268:     ! 
269:     DOUBLE PRECISION, DIMENSION(3) :: CROSS 
270:     DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: A,B 
271:     ! 
272:     CROSS(1) = A(2) * B(3) - A(3) * B(2) 
273:     CROSS(2) = A(3) * B(1) - A(1) * B(3) 
274:     CROSS(3) = A(1) * B(2) - A(2) * B(1) 
275:     ! 
276:   END FUNCTION CROSS 
277:   ! 
278: END MODULE MD 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0