hdiff output

r31974/commons.f90 2017-02-26 20:30:16.266893099 +0000 r31973/commons.f90 2017-02-26 20:30:18.042916621 +0000
 69:      &                 PYSIGNOT, PYEPSNOT, PYA1(3), PYA2(3), PYDPMU, PYDPEPS, PYDPFCT, PYGRAVITYC1, PYGRAVITYC2, & 69:      &                 PYSIGNOT, PYEPSNOT, PYA1(3), PYA2(3), PYDPMU, PYDPEPS, PYDPFCT, PYGRAVITYC1, PYGRAVITYC2, &
 70:      &                 LWRCUT, LWCNSTA, LWCNSTB, LWRCUTSQ, LWRCUT2SQ, DELRC, PAPALP, PAPS, PAPCD, PAPEPS, PAPANG1, PAPANG2, & 70:      &                 LWRCUT, LWCNSTA, LWCNSTB, LWRCUTSQ, LWRCUT2SQ, DELRC, PAPALP, PAPS, PAPCD, PAPEPS, PAPANG1, PAPANG2, &
 71:      &                 DBEPSBB, DBEPSAB, DBSIGBB, DBSIGAB, DBPMU, EFIELD, YKAPPA, YEPS, GEMRC, MREQ, HSEFF, BEPS, & 71:      &                 DBEPSBB, DBEPSAB, DBSIGBB, DBSIGAB, DBPMU, EFIELD, YKAPPA, YEPS, GEMRC, MREQ, HSEFF, BEPS, &
 72:      &                 EPS11, EPS22, EPS12, MRHO11, MRHO22, MRHO12, REQ11, REQ22, REQ12, & 72:      &                 EPS11, EPS22, EPS12, MRHO11, MRHO22, MRHO12, REQ11, REQ22, REQ12, &
 73:      &                 RHOCC0, RHOCC10, RHOCC20,  RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, & 73:      &                 RHOCC0, RHOCC10, RHOCC20,  RHOHH0, RHOHH10, RHOHH20, RHOCH0, RHOC10H, &
 74:      &                 RHOCH10, RHOC20H, RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ, SYMFCTR, & 74:      &                 RHOCH10, RHOC20H, RHOCH20, ALPHACC, ALPHAHH, ALPHACH, DC6CC, DC6HH, DC6CH, KKJ, CCKJ, SYMFCTR, &
 75:      &                 GBANISOTROPYR,GBWELLDEPTHR,PARAMa1,PARAMb1,PARAMc1,PARAMa2,PARAMB2,PARAMc2,PSIGMA0(2),PEPSILON0,& 75:      &                 GBANISOTROPYR,GBWELLDEPTHR,PARAMa1,PARAMb1,PARAMc1,PARAMa2,PARAMB2,PARAMc2,PSIGMA0(2),PEPSILON0,&
 76:      &                 PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2),PYA11(3),PYA21(3),PYA12(3),PYA22(3), & 76:      &                 PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2),PYA11(3),PYA21(3),PYA12(3),PYA22(3), &
 77:      &                 PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, PYCFTHRESH, LJSITECOORDS(3), & 77:      &                 PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, PYCFTHRESH, LJSITECOORDS(3), &
 78:      &                 MSTART,MFINISH,MBSTART1,MBFINISH1,MBSTART2,MBFINISH2,MBHEIGHT1,MBHEIGHT2,ME1,ME2,ME3, & 78:      &                 MSTART,MFINISH,MBSTART1,MBFINISH1,MBSTART2,MBFINISH2,MBHEIGHT1,MBHEIGHT2,ME1,ME2,ME3, &
 79:      &                 BSPTQMAX, BSPTQMIN, PFORCE, CSMNORM, CSMGUIDENORM, CSMEPS, PERCCUT, PERCGROUPCUT, & 79:      &                 BSPTQMAX, BSPTQMIN, PFORCE, CSMNORM, CSMGUIDENORM, CSMEPS, PERCCUT, &
 80:      &                 LOWESTE, PERTSTEP, GCPLUS, & 80:      &                 LOWESTE, PERTSTEP, GCPLUS, &
 81:      &                 KINT, INTFREEZETOL, IMSEPMIN, IMSEPMAX, CONCUTABS, CONCUTFRAC, & 81:      &                 KINT, INTFREEZETOL, IMSEPMIN, IMSEPMAX, CONCUTABS, CONCUTFRAC, &
 82:      &                 LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, & 82:      &                 LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, &
 83:      &                 INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, & 83:      &                 INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, &
 84:      &                 INTLJTOL, INTLJDEL, INTLJEPS, REPCON, INTDGUESS, CHECKREPCUTOFF, INTMINFAC, FREEZETOL, & 84:      &                 INTLJTOL, INTLJDEL, INTLJEPS, REPCON, INTDGUESS, CHECKREPCUTOFF, INTMINFAC, FREEZETOL, &
 85:      &                 LOCALPERMCUT, LOCALPERMCUT2, INTCONCUT, QCIRADSHIFT, MLPLAMBDA, & 85:      &                 LOCALPERMCUT, LOCALPERMCUT2, INTCONCUT, QCIRADSHIFT, MLPLAMBDA, &
 86:      &                 CAPSIDRHO,CAPSIDEPS,SIGMAPENT,RADPENT,SIGMAHEX,RADHEX,SIGMAPH, KLIM, SCA, & 86:      &                 CAPSIDRHO,CAPSIDEPS,SIGMAPENT,RADPENT,SIGMAHEX,RADHEX,SIGMAPH, KLIM, SCA, &
 87:      &                 QCIADDREPCUT, QCIADDREPEPS, MLQLAMBDA 87:      &                 QCIADDREPCUT, QCIADDREPEPS, MLQLAMBDA
 88:  88: 
 89:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, & 89:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, &
 99:      &        ZETT2, GBH_RESTART, RESTART, CONJG, NEWRESTART, AVOID, NATBT, DIFFRACTT, CHRMMT, INTMINT, LB2T, & 99:      &        ZETT2, GBH_RESTART, RESTART, CONJG, NEWRESTART, AVOID, NATBT, DIFFRACTT, CHRMMT, INTMINT, LB2T, &
100:      &        PTMC, BINSTRUCTURES, PROGRESS, MODEL1T, NEWRESTART_MD, CHANGE_TEMP, NOCISTRANS, CHECKCHIRALITY, &100:      &        PTMC, BINSTRUCTURES, PROGRESS, MODEL1T, NEWRESTART_MD, CHANGE_TEMP, NOCISTRANS, CHECKCHIRALITY, &
101:      &        GBT, GBDT, GBDPT, GEMT, LINRODT, RADIFT, CAPBINT, DBPT, DBPTDT, DMBLMT, DMBLPYT, EFIELDT, PAHAT, STOCKAAT, MORSEDPT, &101:      &        GBT, GBDT, GBDPT, GEMT, LINRODT, RADIFT, CAPBINT, DBPT, DBPTDT, DMBLMT, DMBLPYT, EFIELDT, PAHAT, STOCKAAT, MORSEDPT, &
102:      &        MSGBT, MSTBINT, MMRSDPT, MSSTOCKT, LWOTPT, NCAPT, NPAHT, NTIPT, PAHW99T, ELLIPSOIDT, GAYBERNET,&102:      &        MSGBT, MSTBINT, MMRSDPT, MSSTOCKT, LWOTPT, NCAPT, NPAHT, NTIPT, PAHW99T, ELLIPSOIDT, GAYBERNET,&
103:      &        MULTPAHAT, PAPT, PAPBINT, PAPJANT, PTSTSTT, SHIFTED, SILANET, TDHDT, DDMT, WATERDCT, WATERKZT, CHECKDT, CHECKMARKOVT,&103:      &        MULTPAHAT, PAPT, PAPBINT, PAPJANT, PTSTSTT, SHIFTED, SILANET, TDHDT, DDMT, WATERDCT, WATERKZT, CHECKDT, CHECKMARKOVT,&
104:      &        TETHER, HISTSMOOTH, VISITPROP, ARMT, FixedEndMoveT, FIXCOM, RESTORET, QUADT, AMHT, MOVESHELLT, QDT, QD2T, &104:      &        TETHER, HISTSMOOTH, VISITPROP, ARMT, FixedEndMoveT, FIXCOM, RESTORET, QUADT, AMHT, MOVESHELLT, QDT, QD2T, &
105:      &        PARAMONOVPBCX, PARAMONOVPBCY, PARAMONOVPBCZ, PARAMONOVCUTOFF, GAYBERNEDCT, UNFREEZERES, FREEZEALL, &105:      &        PARAMONOVPBCX, PARAMONOVPBCY, PARAMONOVPBCZ, PARAMONOVCUTOFF, GAYBERNEDCT, UNFREEZERES, FREEZEALL, &
106:      &        PROJIT, PROJIHT, LEAPDIHE, DUMPQUT, DUMPBESTT, LJSITE, BLJSITE, LJSITEATTR, DUMPSTEPST, &106:      &        PROJIT, PROJIHT, LEAPDIHE, DUMPQUT, DUMPBESTT, LJSITE, BLJSITE, LJSITEATTR, DUMPSTEPST, &
107:      &        AMBERT, RANDOMSEEDT, PYGPERIODICT, LJCAPSIDT, PYBINARYT, SWAPMOVEST, MOVABLEATOMST, LIGMOVET,DUMPSTRUCTURES, &107:      &        AMBERT, RANDOMSEEDT, PYGPERIODICT, LJCAPSIDT, PYBINARYT, SWAPMOVEST, MOVABLEATOMST, LIGMOVET,DUMPSTRUCTURES, &
108:      &        LJSITECOORDST, VGW, ACKLANDT, G46, DF1T, PULLT, LOCALSAMPLET, CSMT, A9INTET, INTERESTORE, COLDFUSION, &108:      &        LJSITECOORDST, VGW, ACKLANDT, G46, DF1T, PULLT, LOCALSAMPLET, CSMT, A9INTET, INTERESTORE, COLDFUSION, &
109:      &        CSMGUIDET, MULTISITEPYT, CHAPERONINT, AVOIDRESEEDT, OHCELLT, UNFREEZEFINALQ, PERCOLATET, PERCT, PERCACCEPTED, PERCCOMPMARKOV, PERCGROUPT, &109:      &        CSMGUIDET, MULTISITEPYT, CHAPERONINT, AVOIDRESEEDT, OHCELLT, UNFREEZEFINALQ, PERCOLATET, PERCT, PERCACCEPTED, PERCCOMPMARKOV, &
110:      &        GENALT, MINDENSITYT, RESTRICTREGION, RESTRICTREGIONTEST, RESTRICTCYL, ACK1, ACK2, HARMONICF,&110:      &        GENALT, MINDENSITYT, RESTRICTREGION, RESTRICTREGIONTEST, RESTRICTCYL, ACK1, ACK2, HARMONICF,&
111:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&111:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&
112:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, POLYT, SANDBOXT, &112:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, POLYT, SANDBOXT, &
113:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &113:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
114:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &114:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
115:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &115:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
116:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &116:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &
117:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, &117:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, &
118:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &118:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &
119:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &119:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
217:       DOUBLE PRECISION, ALLOCATABLE :: HBONDBESTCOORDS(:,:)217:       DOUBLE PRECISION, ALLOCATABLE :: HBONDBESTCOORDS(:,:)
218:       DOUBLE PRECISION, ALLOCATABLE :: HBONDMARKOV(:), HBONDBEST(:),HBONDQUE(:)218:       DOUBLE PRECISION, ALLOCATABLE :: HBONDMARKOV(:), HBONDBEST(:),HBONDQUE(:)
219:       DOUBLE PRECISION, ALLOCATABLE :: HBONDMAXE(:)219:       DOUBLE PRECISION, ALLOCATABLE :: HBONDMAXE(:)
220:       INTEGER :: HBONDNRES, HBONDMATCHN, MAXHBONDGROUPS220:       INTEGER :: HBONDNRES, HBONDMATCHN, MAXHBONDGROUPS
221: ! END OF HBONDMATRIX VARIABLES221: ! END OF HBONDMATRIX VARIABLES
222: 222: 
223: ! sf344> Modular global optimization variables223: ! sf344> Modular global optimization variables
224:       INTEGER :: MODULARNMIN, MODULARNMAX, MODULARNINCR, MODULARLOWEST,MODULARCURRENTN224:       INTEGER :: MODULARNMIN, MODULARNMAX, MODULARNINCR, MODULARLOWEST,MODULARCURRENTN
225:       LOGICAL :: MODULART,MODULARRESTARTT225:       LOGICAL :: MODULART,MODULARRESTARTT
226:       LOGICAL, ALLOCATABLE :: EXCLUDED(:)226:       LOGICAL, ALLOCATABLE :: EXCLUDED(:)
227:       INTEGER, ALLOCATABLE :: atomingroup(:)227:       CHARACTER(LEN=80) :: MODULARRESTARTFILE
228:       CHARACTER(LEN=80) :: MODULARRESTARTFILE  
229: ! sf344> end of modular global optimization variables228: ! sf344> end of modular global optimization variables
230: 229: 
231:       INTEGER  NGBSITE, NPYSITE, PYBINARYTYPE1, PYSWAP(3), MAXINTERACTIONS230:       INTEGER  NGBSITE, NPYSITE, PYBINARYTYPE1, PYSWAP(3), MAXINTERACTIONS
232: 231: 
233:       DOUBLE PRECISION :: DZP1, DZP2, DZP3, DZP4, DZP5, DZP6, DZP7232:       DOUBLE PRECISION :: DZP1, DZP2, DZP3, DZP4, DZP5, DZP6, DZP7
234:       LOGICAL :: FIELDT, OHT, IHT, TDT, D5HT233:       LOGICAL :: FIELDT, OHT, IHT, TDT, D5HT
235:       DOUBLE PRECISION :: FOH, FIH, FTD, FD5H234:       DOUBLE PRECISION :: FOH, FIH, FTD, FD5H
236: 235: 
237:       !ds656> Box centroid to facilitate global236:       !ds656> Box centroid to facilitate global
238:       !       optimisation of supported clusters.237:       !       optimisation of supported clusters.


r31974/dihedralrotation.f90 2017-02-26 20:30:16.482895959 +0000 r31973/dihedralrotation.f90 2017-02-26 20:30:18.262919534 +0000
128:                   write(myunit,'(A,A,A,F20.10,A,F20.10,A,2F20.10,F10.1)') 'dihedralrotation> rotating group ',trim(adjustl(dihedralgroupnames(i1))),' by ',degreeangle,' from '&128:                   write(myunit,'(A,A,A,F20.10,A,F20.10,A,2F20.10,F10.1)') 'dihedralrotation> rotating group ',trim(adjustl(dihedralgroupnames(i1))),' by ',degreeangle,' from '&
129:      &                  ,angle0*180/pi," new energy ", POTEL,POTELold,IMCSTEP129:      &                  ,angle0*180/pi," new energy ", POTEL,POTELold,IMCSTEP
130:                ENDIF130:                ENDIF
131:             endif131:             endif
132:          enddo132:          enddo
133:       endif133:       endif
134: 134: 
135: 135: 
136: end subroutine dihedralrotstep 136: end subroutine dihedralrotstep 
137: 137: 
138: subroutine dihedralrotation(angle,a0,a1,a2,a3,atomingroup1,stepcoords)138: subroutine dihedralrotation(angle,a0,a1,a2,a3,atomingroup,stepcoords)
139:    USE commons139:    USE commons
140:    IMPLICIT NONE140:    IMPLICIT NONE
141:    INTEGER :: A1, A2, A0,A3,I1141:    INTEGER :: A1, A2, A0,A3,I1
142:    DOUBLE PRECISION :: BVECTOR(3), LENGTH, ANGLE, DUMMYMAT(3,3)=0.0D0, ROTMAT(3,3)142:    DOUBLE PRECISION :: BVECTOR(3), LENGTH, ANGLE, DUMMYMAT(3,3)=0.0D0, ROTMAT(3,3)
143:    DOUBLE PRECISION :: GROUPATOM(3), GROUPATOMROT(3),STEPCOORDS(3*NATOMS),v1(3),v2(3),v3(3),alpha143:    DOUBLE PRECISION :: GROUPATOM(3), GROUPATOMROT(3),STEPCOORDS(3*NATOMS),v1(3),v2(3),v3(3),alpha
144:    LOGICAL :: ATOMINGROUP1(NATOMS)144:    LOGICAL :: ATOMINGROUP(NATOMS)
145: 145: 
146:     !WRITE(*,*) "DIHEDRALROTATION> dihedral angle value to rot", angle*180/3.1415146:     !WRITE(*,*) "DIHEDRALROTATION> dihedral angle value to rot", angle*180/3.1415
147: ! STEP 1147: ! STEP 1
148: ! Produce notmalised bond vector corresponding to the rotation axis148: ! Produce notmalised bond vector corresponding to the rotation axis
149: ! A1 and A2 are the atoms defining this vector149: ! A1 and A2 are the atoms defining this vector
150:    BVECTOR(1)=STEPCOORDS(3*A2-2)-STEPCOORDS(3*A1-2)150:    BVECTOR(1)=STEPCOORDS(3*A2-2)-STEPCOORDS(3*A1-2)
151:    BVECTOR(2)=STEPCOORDS(3*A2-1)-STEPCOORDS(3*A1-1)151:    BVECTOR(2)=STEPCOORDS(3*A2-1)-STEPCOORDS(3*A1-1)
152:    BVECTOR(3)=STEPCOORDS(3*A2  )-STEPCOORDS(3*A1  )152:    BVECTOR(3)=STEPCOORDS(3*A2  )-STEPCOORDS(3*A1  )
153: ! Find length   153: ! Find length   
154:    LENGTH=DSQRT(BVECTOR(1)**2 + BVECTOR(2)**2 + BVECTOR(3)**2)154:    LENGTH=DSQRT(BVECTOR(1)**2 + BVECTOR(2)**2 + BVECTOR(3)**2)
166: ! Interface:166: ! Interface:
167: ! SUBROUTINE RMDRVT(P, RM, DRM1, DRM2, DRM3, GTEST)167: ! SUBROUTINE RMDRVT(P, RM, DRM1, DRM2, DRM3, GTEST)
168: ! P is an un-normalised vector you wish to rotate around. Its length equals the desired rotation in radians168: ! P is an un-normalised vector you wish to rotate around. Its length equals the desired rotation in radians
169: ! RM will return the 3x3 rotation matrix169: ! RM will return the 3x3 rotation matrix
170: ! DRM1-3 are derivative matricies, not needed here170: ! DRM1-3 are derivative matricies, not needed here
171: ! GTEST is also not needed so set to .FALSE.171: ! GTEST is also not needed so set to .FALSE.
172:    CALL RMDRVT(BVECTOR,ROTMAT,DUMMYMAT,DUMMYMAT,DUMMYMAT,.FALSE.)172:    CALL RMDRVT(BVECTOR,ROTMAT,DUMMYMAT,DUMMYMAT,DUMMYMAT,.FALSE.)
173: ! STEP 4173: ! STEP 4
174: ! Rotate group, one atom at a time. First, translate atom so the pivot (end of bond closest to atom) is at the origin  174: ! Rotate group, one atom at a time. First, translate atom so the pivot (end of bond closest to atom) is at the origin  
175:    DO I1=1,NATOMS175:    DO I1=1,NATOMS
176:       IF (ATOMINGROUP1(I1)) THEN176:       IF (ATOMINGROUP(I1)) THEN
177:          GROUPATOM(1)=STEPCOORDS(3*I1-2)-STEPCOORDS(3*A1-2)177:          GROUPATOM(1)=STEPCOORDS(3*I1-2)-STEPCOORDS(3*A1-2)
178:          GROUPATOM(2)=STEPCOORDS(3*I1-1)-STEPCOORDS(3*A1-1)178:          GROUPATOM(2)=STEPCOORDS(3*I1-1)-STEPCOORDS(3*A1-1)
179:          GROUPATOM(3)=STEPCOORDS(3*I1  )-STEPCOORDS(3*A1  )179:          GROUPATOM(3)=STEPCOORDS(3*I1  )-STEPCOORDS(3*A1  )
180: ! Apply the rotation matrix180: ! Apply the rotation matrix
181:          GROUPATOMROT=MATMUL(ROTMAT,GROUPATOM)181:          GROUPATOMROT=MATMUL(ROTMAT,GROUPATOM)
182: ! Translate back to the origin and copy to COORDS182: ! Translate back to the origin and copy to COORDS
183:          STEPCOORDS(3*I1-2)=GROUPATOMROT(1)+STEPCOORDS(3*A1-2)183:          STEPCOORDS(3*I1-2)=GROUPATOMROT(1)+STEPCOORDS(3*A1-2)
184:          STEPCOORDS(3*I1-1)=GROUPATOMROT(2)+STEPCOORDS(3*A1-1)184:          STEPCOORDS(3*I1-1)=GROUPATOMROT(2)+STEPCOORDS(3*A1-1)
185:          STEPCOORDS(3*I1  )=GROUPATOMROT(3)+STEPCOORDS(3*A1  )185:          STEPCOORDS(3*I1  )=GROUPATOMROT(3)+STEPCOORDS(3*A1  )
186:       ENDIF186:       ENDIF


r31974/gay-berne.f90 2017-02-26 20:30:16.734899298 +0000 r31973/gay-berne.f90 2017-02-26 20:30:18.482922450 +0000
4094: 4094: 
4095: 4095: 
4096: !     ----------------------------------------------------------------------------------------------4096: !     ----------------------------------------------------------------------------------------------
4097: 4097: 
4098:       SUBROUTINE TAKESTEPELPSD (NP)4098:       SUBROUTINE TAKESTEPELPSD (NP)
4099: 4099: 
4100: !     THIS ROUTINE TAKES STEP FOR SINGLE-SITE ELLIPSOIDAL BODIES ENSURING NO OVERLAP4100: !     THIS ROUTINE TAKES STEP FOR SINGLE-SITE ELLIPSOIDAL BODIES ENSURING NO OVERLAP
4101: 4101: 
4102:       USE COMMONS 4102:       USE COMMONS 
4103:       use pymodule, only : angle,angle2,twopi4103:       use pymodule, only : angle,angle2,twopi
4104:       use modamber9, only : ligtransstep, ligrotscale 
4105:       use rotations 
4106: 4104: 
4107:       IMPLICIT NONE4105:       IMPLICIT NONE
4108: 4106: 
4109:       INTEGER          :: NP, JMAX, JMAX2, REALNATOMS, OFFSET4107:       INTEGER          :: NP, JMAX, JMAX2, REALNATOMS, OFFSET
4110:       INTEGER          :: J1, J2, J3, J4, J5, J64108:       INTEGER          :: J1, J2, J3, J4, J5, J6
4111: !      INTEGER          :: PTINDX, I, J, K4109: !      INTEGER          :: PTINDX, I, J, K
4112:       LOGICAL          :: OVRLPT4110:       LOGICAL          :: OVRLPT
4113: 4111: 
4114:       DOUBLE PRECISION :: PI, DUMMY, DUMMY2, SAVECOORDS(3*NATOMS), Y(3*NATOMS)4112:       DOUBLE PRECISION :: PI, DUMMY, DUMMY2, SAVECOORDS(3*NATOMS)
4115:       DOUBLE PRECISION :: DIST(3*NATOMS/2), XMASS, YMASS, ZMASS, DMAX, VMAX, VMAX24113:       DOUBLE PRECISION :: DIST(3*NATOMS/2), XMASS, YMASS, ZMASS, DMAX, VMAX, VMAX2
4116:       DOUBLE PRECISION :: VMIN, CMMAX, CMDIST(NATOMS/2), LOCALSTEP4114:       DOUBLE PRECISION :: VMIN, CMMAX, CMDIST(NATOMS/2), LOCALSTEP
4117:       DOUBLE PRECISION :: DPRAND, RANDOM, THETA, PHI4115:       DOUBLE PRECISION :: DPRAND, RANDOM, THETA, PHI
4118:       DOUBLE PRECISION :: AE(3,3), BE(3,3) 4116:       DOUBLE PRECISION :: AE(3,3), BE(3,3) 
4119:       DOUBLE PRECISION :: RI(3), RJ(3), RIJ(3)4117:       DOUBLE PRECISION :: RI(3), RJ(3), RIJ(3)
4120:       DOUBLE PRECISION :: P1(3), P2(3), ABSRIJ, RCUT, XMIN, FMIN, ECFVAL, ligandcoords(3), ligandcoordsrotated(3)4118:       DOUBLE PRECISION :: P1(3), P2(3), ABSRIJ, RCUT, XMIN, FMIN, ECFVAL
4121: ! sf344 additions4119: ! sf344 additions
4122:       LOGICAL          :: DISSOCIATED(NATOMS/2),DISSOC4120:       LOGICAL          :: DISSOCIATED(NATOMS/2),DISSOC
4123:       INTEGER          :: CLOSESTATOMINDEX, K1, K2, i,j, k4121:       INTEGER          :: CLOSESTATOMINDEX, K1, K2
4124:       DOUBLE PRECISION :: MINDISTANCE, x1, x2, y1, y2, z1, z2, RCUTARRAY(NATOMS/2), RVEC(NATOMS/2,NATOMS/2)4122:       DOUBLE PRECISION :: MINDISTANCE, x1, x2, y1, y2, z1, z2, RCUTARRAY(NATOMS/2), RVEC(NATOMS/2,NATOMS/2)
4125:       DOUBLE PRECISION, ALLOCATABLE :: atomblockcentre(:,:), rotationmatrix(:,:,:)4123: !      DOUBLE PRECISION, ALLOCATABLE :: RVEC(:)
4126:       double precision :: st,ct,sph,cph,sps,cosps, randomx, randomy, randomz, randomphi, randomtheta, randompsi 
4127: 4124: 
4128: !      WRITE(*,*) 'NATOMS in takestepelpsd ', NATOMS4125: !      WRITE(*,*) 'NATOMS in takestepelpsd ', NATOMS
4129: !      IF(ALLOCATED(RVEC)) DEALLOCATE(RVEC)4126: !      IF(ALLOCATED(RVEC)) DEALLOCATE(RVEC)
 4127: 
4130:       PI         = 4.D0*ATAN(1.0D0)4128:       PI         = 4.D0*ATAN(1.0D0)
4131:       IF (GBT .OR. GBDT) THEN4129:       IF (GBT .OR. GBDT) THEN
4132: 4130: 
4133:          RCUT       = GBEPSNOT*MAX(GBKAPPA,1.D0)4131:          RCUT       = GBEPSNOT*MAX(GBKAPPA,1.D0)
4134: 4132: 
4135:       ELSE IF (PYGPERIODICT.OR.PYBINARYT) THEN4133:       ELSE IF (PYGPERIODICT.OR.PYBINARYT) THEN
4136:         DO J1=1,NATOMS/24134:         DO J1=1,NATOMS/2
4137:          RCUTARRAY(J1) = 2.0D0 * MAXVAL(PYA1bin(J1,:)) 4135:          RCUTARRAY(J1) = 2.0D0 * MAXVAL(PYA1bin(J1,:)) 
4138:         END DO4136:         END DO
4139:          RCUT = MAXVAL(RCUTARRAY)4137:          RCUT = MAXVAL(RCUTARRAY)
4140: !        WRITE(*,*) 'RCUT=', RCUT 
4141:       ELSE IF (LJCAPSIDT) THEN4138:       ELSE IF (LJCAPSIDT) THEN
4142:         PYA1bin(:,:)=0.5D04139:         PYA1bin(:,:)=0.5D0
4143:         DO J1=1,NATOMS/24140:         DO J1=1,NATOMS/2
4144:          RCUTARRAY(J1) = 2.0D04141:          RCUTARRAY(J1) = 2.0D0
4145:         END DO4142:         END DO
4146:          RCUT = 2.0D04143:          RCUT = 2.0D0
4147:       ENDIF4144:       ENDIF
4148: 4145: 
4149:       REALNATOMS = NATOMS/24146:       REALNATOMS = NATOMS/2
4150:       OFFSET     = 3 * REALNATOMS4147:       OFFSET     = 3 * REALNATOMS
4151: 4148: 
4152:                   SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,NP)4149:                   SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,NP)
4153:         y(:)=COORDS(:,NP) 
4154: ! do the group moves if PERCGROUPT=.TRUE. 
4155:    IF(PERCGROUPT) THEN 
4156:         nblocks=atomingroup(REALNATOMS) !the last number in the atomingroup array is the total number of groups as well 
4157:         WRITE(MYUNIT,*) 'takestepelpsd> taking cooperative rotational and translational steps for ', nblocks, ' blocks'  
4158:         if(allocated(rotationmatrix)) deallocate(rotationmatrix) 
4159:         if(allocated(atomblockcentre)) deallocate(atomblockcentre) 
4160:         allocate(atomblockcentre(nblocks,3),rotationmatrix(nblocks,3,3)) 
4161: !WRITE(MYUNIT,*) 'before rotation and translation:' 
4162: !WRITE(MYUNIT,*) y(:) 
4163:  
4164:  
4165: ! random rotation matrices for each group 
4166:      IF(ligrotscale>0.0D0) then 
4167:        DO i=1,nblocks   
4168:         randomphi=(DPRAND()-0.5)*twopi*ligrotscale 
4169:         randomtheta=(DPRAND()-0.5)*pi*ligrotscale 
4170:         randompsi=(DPRAND()-0.5)*twopi*ligrotscale 
4171:         st=sin(randomtheta) 
4172:         ct=cos(randomtheta) 
4173:         sph=sin(randomphi) 
4174:         cph=cos(randomphi) 
4175:         sps=sin(randompsi) 
4176:         cosps=cos(randompsi) 
4177:  
4178:  
4179:         rotationmatrix(i,1,1)=cosps*cph-ct*sph*sps 
4180:         rotationmatrix(i,2,1)=cosps*sph+ct*cph*sps 
4181:         rotationmatrix(i,3,1)=sps*st 
4182:         rotationmatrix(i,1,2)=-sps*cph-ct*sph*cosps 
4183:         rotationmatrix(i,2,2)=-sps*sph+ct*cph*cosps 
4184:         rotationmatrix(i,3,2)=cosps*st 
4185:         rotationmatrix(i,1,3)=st*sph 
4186:         rotationmatrix(i,2,3)=-st*cph 
4187:         rotationmatrix(i,3,3)=ct 
4188:        END DO 
4189:  
4190:           atomblockcentre(:,:)=0.0D0 
4191:           do i=1,nblocks 
4192:              j=0 
4193:              do k=1,REALNATOMS 
4194:                if(i==atomingroup(k)) then 
4195:                 j=j+1 ! counter of atoms in each block 
4196:                 atomblockcentre(i,1)=atomblockcentre(i,1)+y(3*k-2) 
4197:                 atomblockcentre(i,2)=atomblockcentre(i,2)+y(3*k-1) 
4198:                 atomblockcentre(i,3)=atomblockcentre(i,3)+y(3*k  ) 
4199:                end if 
4200:              end do 
4201: !                write(MYUNIT,*) 'rotation matrix:', rotationmatrix(i,:,:) 
4202:                 atomblockcentre(i,:)=atomblockcentre(i,:)/j 
4203:           end do 
4204:           do i=1,nblocks 
4205:              do k=1,REALNATOMS 
4206:                if(i==atomingroup(k)) then 
4207: !                write(*,*) 'i, k=', i, k 
4208: ! move each block to the origin 
4209:                 ligandcoords(1)=y(3*k-2)-atomblockcentre(i,1) 
4210:                 ligandcoords(2)=y(3*k-1)-atomblockcentre(i,2) 
4211:                 ligandcoords(3)=y(3*k  )-atomblockcentre(i,3) 
4212: ! rotate the block of atoms with random rotation matrix 
4213:                 ligandcoordsrotated=MATMUL(rotationmatrix(i,:,:),ligandcoords) 
4214: !                WRITE(MYUNIT,*) 'ligandcoordsrotated: ', i, k, ligandcoordsrotated(:) 
4215: !                WRITE(MYUNIT,*) 'atomblockcentre:', atomblockcentre(i,:) 
4216: ! translate back the block of atoms to the original centre of their coordinates 
4217:                 y(3*k-2)=ligandcoordsrotated(1)+atomblockcentre(i,1) 
4218:                 y(3*k-1)=ligandcoordsrotated(2)+atomblockcentre(i,2) 
4219:                 y(3*k  )=ligandcoordsrotated(3)+atomblockcentre(i,3) 
4220:                 y(3*k+3*REALNATOMS-2:3*k+3*REALNATOMS)=rot_rotate_aa(y(3*k+3*REALNATOMS-2:3*k+3*REALNATOMS),rot_mx2aa(rotationmatrix(i,:,:))) 
4221:                end if 
4222:              end do 
4223:           end do 
4224: !                WRITE(MYUNIT,*) 'after rotation:' 
4225: !                WRITE(MYUNIT,*) y(:) 
4226:      END IF !ligrotscale 
4227: ! translate each group 
4228:         LOCALSTEP=1.0D0 
4229:            k=1 
4230:            do i=1,nblocks 
4231:              randomx=DPRAND() 
4232:              randomy=DPRAND() 
4233:              randomz=DPRAND() 
4234:              do k=1,REALNATOMS 
4235:               if(i==atomingroup(k)) then 
4236:                  y(3*k-2)=y(3*k-2)+2*randomx*ligtransstep*LOCALSTEP-ligtransstep*LOCALSTEP 
4237:                  y(3*k-1)=y(3*k-1)+2*randomy*ligtransstep*LOCALSTEP-ligtransstep*LOCALSTEP 
4238:                  y(3*k  )=y(3*k  )+2*randomz*ligtransstep*LOCALSTEP-ligtransstep*LOCALSTEP 
4239:               end if 
4240:              end do 
4241:            end do 
4242:  
4243:  
4244:       END IF ! PERCGROUPT 
4245:       COORDS(:,NP)=y(:) 
4246: !WRITE(MYUNIT,*) 'rotated and translated coordinates:' 
4247: !WRITE(MYUNIT,*) y(:) 
4248: 4150: 
4249:       DO J1 = 1,REALNATOMS4151:       DO J1 = 1,REALNATOMS
4250: 4152: 
4251:          J2     = 3*J14153:          J2     = 3*J1
4252:         IF (PARAMONOVPBCX) THEN4154:         IF (PARAMONOVPBCX) THEN
4253: !ensure x component of particle 1 vector is within BoxLx/2 of zero.4155: !ensure x component of particle 1 vector is within BoxLx/2 of zero.
4254: !if it isn't then subtract integer number of boxlx's such that it is.4156: !if it isn't then subtract integer number of boxlx's such that it is.
4255:                 COORDS(J2-2,NP)=COORDS(J2-2,NP)-BOXLX*NINT(COORDS(J2-2,NP)/BOXLX)4157:                 COORDS(J2-2,NP)=COORDS(J2-2,NP)-BOXLX*NINT(COORDS(J2-2,NP)/BOXLX)
4256:         ENDIF4158:         ENDIF
4257: 4159: 
4446:                         4348:                         
4447:                IF (ABSRIJ < RCUT) THEN 4349:                IF (ABSRIJ < RCUT) THEN 
4448: 4350: 
4449: !     DETERMINE ELLIPTIC CONTACT FUNCTION4351: !     DETERMINE ELLIPTIC CONTACT FUNCTION
4450: 4352: 
4451:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE, BE, RIJ, XMIN, FMIN)4353:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE, BE, RIJ, XMIN, FMIN)
4452: 4354: 
4453:                   ECFVAL = - FMIN4355:                   ECFVAL = - FMIN
4454: ! allow for a slight overlap 4356: ! allow for a slight overlap 
4455:                   IF (ECFVAL < PYOVERLAPTHRESH) THEN4357:                   IF (ECFVAL < PYOVERLAPTHRESH) THEN
4456:                         WRITE(MYUNIT,*) 'atoms overlapping', J1, J2, ECFVAL, ABSRIJ4358: !                        WRITE(*,*) 'atoms overlapping', J1, J2, ECFVAL, ABSRIJ
4457:                      OVRLPT = .TRUE.4359:                      OVRLPT = .TRUE.
4458:                      GO TO 954360:                      GO TO 95
4459:                   ENDIF4361:                   ENDIF
4460: 4362: 
4461:                ENDIF4363:                ENDIF
4462: 4364: 
4463:             ENDDO  ! END LOOP OVER J24365:             ENDDO  ! END LOOP OVER J2
4464:          ENDDO  ! END WHILE4366:          ENDDO  ! END WHILE
4465:       ENDDO  ! END LOOP OVER J1   4367:       ENDDO  ! END LOOP OVER J1   
4466: 4368: 
4467:       IF (FREEZE) THEN4369:                      IF (FREEZE) THEN
4468:          DO J2=1,NATOMS4370:                         DO J2=1,NATOMS
4469:            IF (FROZEN(J2)) THEN4371:                            IF (FROZEN(J2)) THEN
4470:               COORDS(3*(J2-1)+1:3*(J2-1)+3,NP)=SAVECOORDS(3*(J2-1)+1:3*(J2-1)+3)4372:                               COORDS(3*(J2-1)+1:3*(J2-1)+3,NP)=SAVECOORDS(3*(J2-1)+1:3*(J2-1)+3)
4471:            ENDIF4373:                            ENDIF
4472:          ENDDO4374:                         ENDDO
4473:       ENDIF4375:                      ENDIF
 4376: 
 4377: 
 4378: ! now determine whether any particle has dissociated
 4379:  IF(PYGPERIODICT) THEN
 4380:   DO
 4381:   DISSOC=.FALSE.
 4382:    DO J1=1,natoms/2-1
 4383:     J2=3*J1
 4384:     x1=COORDS(J2-2,NP)
 4385:     y1=COORDS(J2-1,NP)
 4386:     z1=COORDS(J2,NP)
 4387:    
 4388:     DISSOCIATED(J1)=.TRUE.
 4389:     DO K1=1,natoms/2
 4390:       K2=3*K1
 4391:       x2=COORDS(K2-2,NP)
 4392:       y2=COORDS(K2-1,NP)
 4393:       z2=COORDS(K2,NP)
 4394:       RVEC(J1,K1)=SQRT((x2-x1)**2+(y2-y1)**2+(z2-z1)**2)
 4395: !      WRITE(*,*) J1,K1,RVEC(J1,K1), 4*RCUT
 4396:       IF(RVEC(J1,K1)<4*RCUT) DISSOCIATED(J1)=.FALSE.
 4397: !      WRITE(*,*) 'in here, cutoff=',pcutoff,J1
 4398:     END DO
 4399:     IF(DISSOCIATED(J1)) THEN
 4400:         ! determine the closest neighbor
 4401:         MINDISTANCE = HUGE(1.0D0)
 4402:         DO K1=1,natoms/2
 4403:          IF(K1==J1) CYCLE
 4404:          IF(RVEC(J1,K1)<MINDISTANCE) THEN
 4405:            MINDISTANCE=RVEC(J1,K1)
 4406:            CLOSESTATOMINDEX=K1
 4407:          END IF
 4408:         END DO
 4409:         WRITE(MYUNIT,'(A,I6,A,F10.6,I6)') 'Atom ', J1, ' is dissociated, distance to the closest atom: ', &
 4410:         & MINDISTANCE, CLOSESTATOMINDEX
 4411:         DISSOC=.TRUE.
 4412:         ! move the dissociated atom closer to the closest atom
 4413:         COORDS(J2-2,NP)=0.5D0*(COORDS(3*CLOSESTATOMINDEX-2,NP)+COORDS(J2-2,NP))
 4414:         COORDS(J2-1,NP)=0.5D0*(COORDS(3*CLOSESTATOMINDEX-1,NP)+COORDS(J2-1,NP))
 4415:         COORDS(J2  ,NP)=0.5D0*(COORDS(3*CLOSESTATOMINDEX  ,NP)+COORDS(J2  ,NP))
 4416:     END IF
 4417:    END DO
 4418:    IF(.NOT.DISSOC) EXIT
 4419:   END DO
 4420:         DO J1=NATOMS/2+1,NATOMS
 4421:             J2=3*J1
 4422:             angle2=dot_product(COORDS(J2-2:J2,NP),COORDS(J2-2:J2,NP))
 4423: 
 4424:             if(angle2>twopi**2) then
 4425:                 angle2=sqrt(angle2)
 4426:                 angle = angle2 - int(angle2/twopi)*twopi
 4427:                 COORDS(J2-2:J2,NP)=COORDS(J2-2:J2,NP)/angle2 * angle
 4428:             end if
 4429:         END DO
4474: 4430: 
4475: RETURN4431: RETURN
 4432: !     CHECK FOR OVERLAP AGAIN (not working!)
 4433:       DO J1 = 1, REALNATOMS
4476: 4434: 
 4435:          J3 = 3*J1
 4436:          J5 = OFFSET + J3
4477: 4437: 
4478: !! now determine whether any particle has dissociated4438:          OVRLPT = .TRUE.
4479: ! IF(PYGPERIODICT) THEN4439: 
4480: !  DO4440: 195       DO WHILE (OVRLPT)
4481: !  DISSOC=.FALSE.4441: 
4482: !   DO J1=1,natoms/2-14442:             LOCALSTEP = 0.0D0
4483: !    J2=3*J14443:             IF (TMOVE(NP)) LOCALSTEP = STEP(NP)
4484: !    x1=COORDS(J2-2,NP)4444: 
4485: !    y1=COORDS(J2-1,NP)4445:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0
4486: !    z1=COORDS(J2,NP)4446:             COORDS(J3-2,NP) = COORDS(J3-2,NP) + LOCALSTEP*RANDOM
4487: !   4447:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0
4488: !    DISSOCIATED(J1)=.TRUE.4448:             COORDS(J3-1,NP) = COORDS(J3-1,NP) + LOCALSTEP*RANDOM
4489: !    DO K1=J1+1,natoms/24449:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0
4490: !      K2=3*K14450:             COORDS(J3,NP)   = COORDS(J3,NP) + LOCALSTEP*RANDOM
4491: !      x2=COORDS(K2-2,NP)4451: 
4492: !      y2=COORDS(K2-1,NP)4452:             LOCALSTEP = 0.0D0
4493: !      z2=COORDS(K2,NP)4453:             IF (OMOVE(NP)) LOCALSTEP = OSTEP(NP)
4494: !      RVEC(J1,K1)=SQRT((x2-x1)**2+(y2-y1)**2+(z2-z1)**2)4454: 
4495: !!      WRITE(*,*) J1,K1,RVEC(J1,K1), 4*RCUT4455:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0
4496: !      IF(RVEC(J1,K1)<4*RCUT) DISSOCIATED(J1)=.FALSE.4456:             COORDS(J5-2,NP) = COORDS(J5-2,NP) + LOCALSTEP*RANDOM
4497: !!      WRITE(*,*) 'in here, cutoff=',pcutoff,J14457:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0
4498: !    END DO4458:             COORDS(J5-1,NP) = COORDS(J5-1,NP) + LOCALSTEP*RANDOM
4499: !    IF(DISSOCIATED(J1)) THEN4459:             RANDOM          = (DPRAND() - 0.5D0)*2.0D0
4500: !        ! determine the closest neighbor4460:             COORDS(J5,NP)   = COORDS(J5,NP) + LOCALSTEP*RANDOM
4501: !        MINDISTANCE = HUGE(1.0D0)4461:           
4502: !        DO K1=1,natoms/24462:             OVRLPT = .FALSE.
4503: !         IF(K1==J1) CYCLE4463: 
4504: !         IF(RVEC(J1,K1)<MINDISTANCE) THEN4464:             RI(:)  = COORDS(J3-2:J3,NP)
4505: !           MINDISTANCE=RVEC(J1,K1)4465:             P1(:)  = COORDS(J5-2:J5,NP)
4506: !           CLOSESTATOMINDEX=K14466: 
4507: !         END IF4467: !     ROTATION MATRIX
4508: !        END DO4468: 
4509: !        WRITE(MYUNIT,'(A,I6,A,F10.6,I6,F10.6)') 'Atom ', J1, ' is dissociated, distance to the closest atom: ', &4469:             CALL ELPSDRTMT (P1, ESA, AE)
4510: !        & MINDISTANCE, CLOSESTATOMINDEX, RCUT4470: 
4511: !        DISSOC=.TRUE.4471:             DO J2 = 1, REALNATOMS
4512: !        ! move the dissociated atom closer to the closest atom4472: 
4513: !        COORDS(J2-2,NP)=0.5D0*(COORDS(3*CLOSESTATOMINDEX-2,NP)+COORDS(J2-2,NP))4473:                IF (J2 == J1) CYCLE
4514: !        COORDS(J2-1,NP)=0.5D0*(COORDS(3*CLOSESTATOMINDEX-1,NP)+COORDS(J2-1,NP))4474:  
4515: !        COORDS(J2  ,NP)=0.5D0*(COORDS(3*CLOSESTATOMINDEX  ,NP)+COORDS(J2  ,NP))4475:                J4    = 3*J2
4516: !    END IF4476:                J6    = OFFSET + J4
4517: !   END DO4477:                RJ(:) = COORDS(J4-2:J4,NP)
4518: !   IF(.NOT.DISSOC) EXIT4478:                P2(:) = COORDS(J6-2:J6,NP)
4519: !  END DO4479: 
4520: !        DO J1=NATOMS/2+1,NATOMS4480: !     ROTATION MATRIX
4521: !            J2=3*J14481: 
4522: !            angle2=dot_product(COORDS(J2-2:J2,NP),COORDS(J2-2:J2,NP))!4482:                CALL ELPSDRTMT (P2, ESA, BE)
4523: !4483: 
4524: !            if(angle2>twopi**2) then4484:                RIJ    = RI - RJ
4525: !                angle2=sqrt(angle2)4485:                ABSRIJ = DSQRT(DOT_PRODUCT(RIJ,RIJ))
4526: !                angle = angle2 - int(angle2/twopi)*twopi4486:                         
4527: !                COORDS(J2-2:J2,NP)=COORDS(J2-2:J2,NP)/angle2 * angle4487:                IF (ABSRIJ < RCUT) THEN 
4528: !            end if4488: 
4529: !        END DO4489: !     DETERMINE ELLIPTIC CONTACT FUNCTION
 4490: 
 4491:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE, BE, RIJ, XMIN, FMIN)
 4492: 
 4493:                   ECFVAL = - FMIN
 4494:  
 4495:                   IF (ECFVAL < 1.D0) THEN
 4496:                      OVRLPT = .TRUE.
 4497:                      GO TO 195
 4498:                   ENDIF
 4499: 
 4500:                ENDIF
 4501: 
 4502:             ENDDO  ! END LOOP OVER J2
 4503: 
 4504:          ENDDO  ! END WHILE
 4505:  
 4506:       ENDDO  ! END LOOP OVER J1   
4530: 4507: 
4531: !RETURN4508:  END IF
4532: 4509: 
4533:       END SUBROUTINE TAKESTEPELPSD4510:       END SUBROUTINE TAKESTEPELPSD
4534: 4511: 
4535: SUBROUTINE TAKESTEPSWAPMOVES(NP)4512: SUBROUTINE TAKESTEPSWAPMOVES(NP)
4536: use commons, only : NATOMS, COORDS, PYBINARYTYPE1,MYUNIT,PYSWAP!,SWAPMOVEST4513: use commons, only : NATOMS, COORDS, PYBINARYTYPE1,MYUNIT,PYSWAP!,SWAPMOVEST
4537: 4514: 
4538: implicit none4515: implicit none
4539: 4516: 
4540: integer :: NP, RANDOM1, RANDOM2, J1, J24517: integer :: NP, RANDOM1, RANDOM2, J1, J2
4541: double precision :: COORDSSTORE(3,4),DPRAND4518: double precision :: COORDSSTORE(3,4),DPRAND


r31974/grouprotation.f90 2017-02-26 20:30:16.950902157 +0000 r31973/grouprotation.f90 2017-02-26 20:30:18.702925361 +0000
130:       ENDIF130:       ENDIF
131: 131: 
132:       END SUBROUTINE GROUPROTSTEP132:       END SUBROUTINE GROUPROTSTEP
133: 133: 
134: ! The GROUPROTATION subroutine allows for almost any rotation of a defined set of atoms.134: ! The GROUPROTATION subroutine allows for almost any rotation of a defined set of atoms.
135: ! The rotation axis is defined by two atoms (BATOMS1 and BATOM2), the group of atoms to 135: ! The rotation axis is defined by two atoms (BATOMS1 and BATOM2), the group of atoms to 
136: ! rotate is defined by the logical array ATOMINGROUP (if element is .TRUE., that atom is136: ! rotate is defined by the logical array ATOMINGROUP (if element is .TRUE., that atom is
137: ! in the group), ANGLE is the rotation angle in radians and STEPCOORDS contains the current137: ! in the group), ANGLE is the rotation angle in radians and STEPCOORDS contains the current
138: ! atomic coordinates.138: ! atomic coordinates.
139: !139: !
140:       SUBROUTINE GROUPROTATION(BATOM1,BATOM2,ANGLE,ATOMINGROUP1,STEPCOORDS)140:       SUBROUTINE GROUPROTATION(BATOM1,BATOM2,ANGLE,ATOMINGROUP,STEPCOORDS)
141:       USE commons141:       USE commons
142:       USE MOVES142:       USE MOVES
143:       IMPLICIT NONE143:       IMPLICIT NONE
144:       INTEGER :: BATOM1, BATOM2, I1144:       INTEGER :: BATOM1, BATOM2, I1
145:       DOUBLE PRECISION :: BVECTOR(3), LENGTH, ANGLE, DUMMYMAT(3,3)=0.0D0, ROTMAT(3,3)145:       DOUBLE PRECISION :: BVECTOR(3), LENGTH, ANGLE, DUMMYMAT(3,3)=0.0D0, ROTMAT(3,3)
146:       DOUBLE PRECISION :: GROUPATOM(3), GROUPATOMROT(3), STEPCOORDS(3*NATOMS)146:       DOUBLE PRECISION :: GROUPATOM(3), GROUPATOMROT(3), STEPCOORDS(3*NATOMS)
147:       LOGICAL :: ATOMINGROUP1(NATOMS)147:       LOGICAL :: ATOMINGROUP(NATOMS)
148: ! ===============================TESTING ROTATE_ABOUT_AXIS============================148: ! ===============================TESTING ROTATE_ABOUT_AXIS============================
149: !      DOUBLE PRECISION, DIMENSION(3*NATOMS)  :: COORDS_COPY149: !      DOUBLE PRECISION, DIMENSION(3*NATOMS)  :: COORDS_COPY
150: !      INTEGER, DIMENSION(:), ALLOCATABLE     :: ATOM_LIST150: !      INTEGER, DIMENSION(:), ALLOCATABLE     :: ATOM_LIST
151: !      INTEGER                                :: NUM_ROTATING_ATOMS151: !      INTEGER                                :: NUM_ROTATING_ATOMS
152: !      INTEGER                                :: I2152: !      INTEGER                                :: I2
153: !      DOUBLE PRECISION, DIMENSION(3)         :: AXIS_START, AXIS_END153: !      DOUBLE PRECISION, DIMENSION(3)         :: AXIS_START, AXIS_END
154: !      DOUBLE PRECISION                       :: PI, ANGLE_DEGREES154: !      DOUBLE PRECISION                       :: PI, ANGLE_DEGREES
155: !155: !
156: !! Take a copy of the coordinates156: !! Take a copy of the coordinates
157: !      COORDS_COPY(:) = STEPCOORDS(:)157: !      COORDS_COPY(:) = STEPCOORDS(:)
213: ! Interface:213: ! Interface:
214: ! SUBROUTINE RMDRVT(P, RM, DRM1, DRM2, DRM3, GTEST)214: ! SUBROUTINE RMDRVT(P, RM, DRM1, DRM2, DRM3, GTEST)
215: ! P is an un-normalised vector you wish to rotate around. Its length equals the desired rotation in radians215: ! P is an un-normalised vector you wish to rotate around. Its length equals the desired rotation in radians
216: ! RM will return the 3x3 rotation matrix216: ! RM will return the 3x3 rotation matrix
217: ! DRM1-3 are derivative matricies, not needed here217: ! DRM1-3 are derivative matricies, not needed here
218: ! GTEST is also not needed so set to .FALSE.218: ! GTEST is also not needed so set to .FALSE.
219:       CALL RMDRVT(BVECTOR,ROTMAT,DUMMYMAT,DUMMYMAT,DUMMYMAT,.FALSE.)219:       CALL RMDRVT(BVECTOR,ROTMAT,DUMMYMAT,DUMMYMAT,DUMMYMAT,.FALSE.)
220: ! STEP 4220: ! STEP 4
221: ! Rotate group, one atom at a time. First, translate atom so the pivot (end of bond closest to atom) is at the origin  221: ! Rotate group, one atom at a time. First, translate atom so the pivot (end of bond closest to atom) is at the origin  
222:       DO I1=1,NATOMS222:       DO I1=1,NATOMS
223:          IF (ATOMINGROUP1(I1)) THEN223:          IF (ATOMINGROUP(I1)) THEN
224:             GROUPATOM(1)=STEPCOORDS(3*I1-2)-STEPCOORDS(3*BATOM2-2)224:             GROUPATOM(1)=STEPCOORDS(3*I1-2)-STEPCOORDS(3*BATOM2-2)
225:             GROUPATOM(2)=STEPCOORDS(3*I1-1)-STEPCOORDS(3*BATOM2-1)225:             GROUPATOM(2)=STEPCOORDS(3*I1-1)-STEPCOORDS(3*BATOM2-1)
226:             GROUPATOM(3)=STEPCOORDS(3*I1  )-STEPCOORDS(3*BATOM2  )226:             GROUPATOM(3)=STEPCOORDS(3*I1  )-STEPCOORDS(3*BATOM2  )
227: ! Apply the rotation matrix227: ! Apply the rotation matrix
228:             GROUPATOMROT=MATMUL(ROTMAT,GROUPATOM)228:             GROUPATOMROT=MATMUL(ROTMAT,GROUPATOM)
229: ! Translate back to the origin and copy to COORDS229: ! Translate back to the origin and copy to COORDS
230:             STEPCOORDS(3*I1-2)=GROUPATOMROT(1)+STEPCOORDS(3*BATOM2-2)230:             STEPCOORDS(3*I1-2)=GROUPATOMROT(1)+STEPCOORDS(3*BATOM2-2)
231:             STEPCOORDS(3*I1-1)=GROUPATOMROT(2)+STEPCOORDS(3*BATOM2-1)231:             STEPCOORDS(3*I1-1)=GROUPATOMROT(2)+STEPCOORDS(3*BATOM2-1)
232:             STEPCOORDS(3*I1  )=GROUPATOMROT(3)+STEPCOORDS(3*BATOM2  )232:             STEPCOORDS(3*I1  )=GROUPATOMROT(3)+STEPCOORDS(3*BATOM2  )
233:          ENDIF233:          ENDIF


r31974/keywords.f 2017-02-26 20:30:17.170905071 +0000 r31973/keywords.f 2017-02-26 20:30:18.922928277 +0000
995:       PULLT=.FALSE.995:       PULLT=.FALSE.
996:       CSMT=.FALSE.996:       CSMT=.FALSE.
997:       CSMGUIDET=.FALSE.997:       CSMGUIDET=.FALSE.
998:       CSMEPS=1.0D-6998:       CSMEPS=1.0D-6
999:       CSMSTEPS=1999:       CSMSTEPS=1
1000:       CSMQUENCHES=11000:       CSMQUENCHES=1
1001:       CSMMAXIT=01001:       CSMMAXIT=0
1002:       CHECKMARKOVT=.FALSE.1002:       CHECKMARKOVT=.FALSE.
1003:       PERCOLATET=.FALSE.1003:       PERCOLATET=.FALSE.
1004:       PERCCUT=1.0D1001004:       PERCCUT=1.0D100
1005:       PERCGROUPCUT=1.0D100 
1006:       PERCGROUPT=.FALSE. 
1007:       GENALT=.FALSE.1005:       GENALT=.FALSE.
1008: !1006: !
1009: ! Reservoir list parameters for PT / BSPT.1007: ! Reservoir list parameters for PT / BSPT.
1010: !1008: !
1011:       RESERVOIRT=.FALSE.1009:       RESERVOIRT=.FALSE.
1012:       USERES=01010:       USERES=0
1013:       EXEQ=01011:       EXEQ=0
1014:       PERTSTEP=0.01D01012:       PERTSTEP=0.01D0
1015:       RES_PSWAP = 0.0D01013:       RES_PSWAP = 0.0D0
1016:       BETA_RES = 0.0D01014:       BETA_RES = 0.0D0
2541:         IF (ligcartstep.gt.0) THEN2539:         IF (ligcartstep.gt.0) THEN
2542:            WRITE(MYUNIT,'(A,G10.3,A)') ' keyword> cartesian perturbations of up to ',ligcartstep,' will be applied to the ligand'2540:            WRITE(MYUNIT,'(A,G10.3,A)') ' keyword> cartesian perturbations of up to ',ligcartstep,' will be applied to the ligand'
2543:         ENDIF2541:         ENDIF
2544: ! sf344> block moves. To be used for rigid body moves of multiple ligands or small clusters of molecules. Rigid body move2542: ! sf344> block moves. To be used for rigid body moves of multiple ligands or small clusters of molecules. Rigid body move
2545: !        parameters will be set using LIGMOVE. Syntax: BLOCKMOVE nblocks block1 block2 ... blockN. The file movableatomslist2543: !        parameters will be set using LIGMOVE. Syntax: BLOCKMOVE nblocks block1 block2 ... blockN. The file movableatomslist
2546: !        has to contain all atom indices, and these will be separated into nblocks parts, each block containing block_i number2544: !        has to contain all atom indices, and these will be separated into nblocks parts, each block containing block_i number
2547: !        of atoms.2545: !        of atoms.
2548:       ELSE IF(WORD.EQ.'BLOCKMOVE') THEN2546:       ELSE IF(WORD.EQ.'BLOCKMOVE') THEN
2549:         BLOCKMOVET=.TRUE.2547:         BLOCKMOVET=.TRUE.
2550:         CALL READI(NBLOCKS)2548:         CALL READI(NBLOCKS)
2551:         IF(ALLOCATED(ATOMSINBLOCK)) DEALLOCATE(ATOMSINBLOCK)2549:         DEALLOCATE(ATOMSINBLOCK)
2552:         ALLOCATE(ATOMSINBLOCK(NBLOCKS))2550:         ALLOCATE(ATOMSINBLOCK(NBLOCKS))
2553:         DO J1=1,NBLOCKS2551:         DO J1=1,NBLOCKS
2554:            CALL READI(ATOMSINBLOCK(J1))2552:            CALL READI(ATOMSINBLOCK(J1))
2555:         END DO2553:         END DO
2556:         WRITE(MYUNIT,*) ' keywords> atoms in movableatomlist will be moved independently in blocks of ', ATOMSINBLOCK(:)2554:         WRITE(MYUNIT,*) ' keywords> atoms in movableatomlist will be moved independently in blocks of ', ATOMSINBLOCK(:)
2557: !2555: !
2558: ! Explicit dump of image coordinates in xyz format for intlbfgs. Should2556: ! Explicit dump of image coordinates in xyz format for intlbfgs. Should
2559: ! be set .TRUE. if DEBUG is set.2557: ! be set .TRUE. if DEBUG is set.
2560: !2558: !
2561:       ELSE IF (WORD == 'DUMPINTXYZ') THEN2559:       ELSE IF (WORD == 'DUMPINTXYZ') THEN
5162:          PBGLUET=.TRUE.5160:          PBGLUET=.TRUE.
5163: 5161: 
5164:       ELSE IF (WORD.EQ.'PERCOLATE') THEN5162:       ELSE IF (WORD.EQ.'PERCOLATE') THEN
5165:          PERCOLATET=.TRUE.5163:          PERCOLATET=.TRUE.
5166:          PERCACCEPTED=.FALSE.5164:          PERCACCEPTED=.FALSE.
5167:          PERCCOMPMARKOV=.FALSE.5165:          PERCCOMPMARKOV=.FALSE.
5168:          IF (NITEMS.GT.1) CALL READF(PERCCUT)5166:          IF (NITEMS.GT.1) CALL READF(PERCCUT)
5169:          IF (NITEMS.GT.2) CALL READF(K_COMP)5167:          IF (NITEMS.GT.2) CALL READF(K_COMP)
5170:          IF (NITEMS.GT.3) CALL READF(GUIDECUT)5168:          IF (NITEMS.GT.3) CALL READF(GUIDECUT)
5171: 5169: 
5172:       ELSE IF (WORD.EQ.'PERCGROUP') THEN 
5173:          PERCGROUPT=.TRUE. 
5174:          CALL READF(PERCGROUPCUT) 
5175:          CALL READF(LIGTRANSSTEP) 
5176:          CALL READF(LIGROTSCALE) 
5177:  
5178:       ELSE IF (WORD.EQ.'PERIODIC') THEN5170:       ELSE IF (WORD.EQ.'PERIODIC') THEN
5179:          PERIODIC=.TRUE.5171:          PERIODIC=.TRUE.
5180:          CALL READF(XX)5172:          CALL READF(XX)
5181:          BOXLX=XX5173:          BOXLX=XX
5182:          BOX3D(1) = XX5174:          BOX3D(1) = XX
5183:          IF (NITEMS.GT.2) THEN5175:          IF (NITEMS.GT.2) THEN
5184:             CALL READF(XX)5176:             CALL READF(XX)
5185:             BOXLY=XX5177:             BOXLY=XX
5186:             BOX3D(2) = XX5178:             BOX3D(2) = XX
5187:             IF (NITEMS.GT.3) THEN5179:             IF (NITEMS.GT.3) THEN


r31974/mc.F 2017-02-26 20:30:17.390907984 +0000 r31973/mc.F 2017-02-26 20:30:19.150931296 +0000
992: !               seed the random number generator with system time + MYNODE (for MPI runs)992: !               seed the random number generator with system time + MYNODE (for MPI runs)
993:                   IF(RANDOMSEEDT) THEN993:                   IF(RANDOMSEEDT) THEN
994:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)994:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)
995:                      itime1= values(6)*60 + values(7)995:                      itime1= values(6)*60 + values(7)
996:                      CALL SDPRND(itime1+MYNODE)996:                      CALL SDPRND(itime1+MYNODE)
997:                   END IF997:                   END IF
998:                 call py_takestep(jp)998:                 call py_takestep(jp)
999:                 if (rigidmdt) call rigidmd_nvt999:                 if (rigidmdt) call rigidmd_nvt
1000: 1000: 
1001:                ELSE IF (PYGPERIODICT.OR.PYBINARYT.OR.LJCAPSIDT) THEN1001:                ELSE IF (PYGPERIODICT.OR.PYBINARYT.OR.LJCAPSIDT) THEN
1002:                 IF(PERCGROUPT) THEN 
1003: ! sf344> group particles according to percolation, reorder them 
1004:                   CALL PERCGROUP(COORDS(1:3*NATOMS,JP),NATOMS,PERCGROUPCUT,DEBUG,MYUNIT,RIGID) 
1005: !                 in this case COORDSO is the same as COORDS, except atoms are being reshuffled. Reset 
1006: !                  COORDSO to get rid of sanity check later on 
1007:                   COORDSO(:,JP)=COORDS(:,JP) 
1008:                 END IF 
1009: ! sf344> particles are now ordered according to the percolated groups they are in. Group indices are 
1010: !        saved in atomingroup(:) 
1011: !               seed the random number generator with system time + MYNODE (for MPI runs)1002: !               seed the random number generator with system time + MYNODE (for MPI runs)
1012:                   IF(RANDOMSEEDT) THEN1003:                   IF(RANDOMSEEDT) THEN
1013:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)1004:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)
1014:                      itime1= values(6)*60 + values(7)1005:                      itime1= values(6)*60 + values(7)
1015:                      CALL SDPRND(itime1+MYNODE)1006:                      CALL SDPRND(itime1+MYNODE)
1016:                   END IF1007:                   END IF
1017:                   IF(SWAPMOVEST) THEN1008:                   IF(SWAPMOVEST) THEN
1018:                      CALL TAKESTEPSWAPMOVES(JP)1009:                      CALL TAKESTEPSWAPMOVES(JP)
1019:                      CALL TAKESTEPELPSD(JP)1010:                      CALL TAKESTEPELPSD(JP)
1020:                   ELSE1011:                   ELSE
1955:                      ELSE1946:                      ELSE
1956:                         EPREV(JP)=HBONDESAVE1947:                         EPREV(JP)=HBONDESAVE
1957: !                        WRITE(MYUNIT,'(A)') " HBONDMATRIX> Markov energy reset to pass sanity checks"1948: !                        WRITE(MYUNIT,'(A)') " HBONDMATRIX> Markov energy reset to pass sanity checks"
1958:                      ENDIF1949:                      ENDIF
1959:                   ENDIF1950:                   ENDIF
1960:                ENDIF1951:                ENDIF
1961: !1952: !
1962: !  Sanity check to make sure the Markov energy agrees with COORDSO. 1953: !  Sanity check to make sure the Markov energy agrees with COORDSO. 
1963: !  Stop if not true.1954: !  Stop if not true.
1964: !1955: !
1965:                IF ((DEBUG.OR.CHECKMARKOVT).AND.(.NOT.(HBONDMATRIX.OR.PERMOPT.OR.PERMINVOPT.OR.FEBHT.OR.QALCST.OR. 1956:                IF ((DEBUG.OR.CHECKMARKOVT).AND.(.NOT.(HBONDMATRIX.OR.PERMOPT.OR.PERMINVOPT.OR.FEBHT.OR.QALCST.OR.MODULART))) THEN
1966: &                               MODULART.OR.PERCGROUPT))) THEN 
1967: ! jwrm2> If we're using percolate, we need to set compression to whether it was on or not when calculating the Markov energy1957: ! jwrm2> If we're using percolate, we need to set compression to whether it was on or not when calculating the Markov energy
1968:                   QUENCHCOMP = COMPON1958:                   QUENCHCOMP = COMPON
1969:                   IF (PERCOLATET) COMPON = PERCCOMPMARKOV1959:                   IF (PERCOLATET) COMPON = PERCCOMPMARKOV
1970:                   CALL POTENTIAL(COORDSO(:,JP),GRAD,OPOTEL,.FALSE.,.FALSE.)1960:                   CALL POTENTIAL(COORDSO(:,JP),GRAD,OPOTEL,.FALSE.,.FALSE.)
1971:                   IF (ABS(OPOTEL-EPREV(JP)).GT.ABS(ECONV)) THEN1961:                   IF (ABS(OPOTEL-EPREV(JP)).GT.ABS(ECONV)) THEN
1972:                      IF (EVAP) THEN1962:                      IF (EVAP) THEN
1973:                         WRITE(MYUNIT,'(3(A,G20.10))') 'mc> WARNING - energy for saved coordinates ',OPOTEL,1963:                         WRITE(MYUNIT,'(3(A,G20.10))') 'mc> WARNING - energy for saved coordinates ',OPOTEL,
1974:      &                     ' differs from Markov energy ',EPREV(JP),' because an atom moved outside the container'1964:      &                     ' differs from Markov energy ',EPREV(JP),' because an atom moved outside the container'
1975:                      ELSE1965:                      ELSE
1976:                         WRITE(MYUNIT,'(2(A,G20.10))') 'mc> ERROR - energy for coordinates in COORDSO=',OPOTEL,1966:                         WRITE(MYUNIT,'(2(A,G20.10))') 'mc> ERROR - energy for coordinates in COORDSO=',OPOTEL,


r31974/perc.f90 2017-02-26 20:30:17.606910846 +0000 r31973/perc.f90 2017-02-26 20:30:19.370934208 +0000
 94:      IF (DEBUG) THEN 94:      IF (DEBUG) THEN
 95:         DO J1=1,NSITES 95:         DO J1=1,NSITES
 96:            WRITE(MYUNIT,'(2I6,3G20.10)') J1,NDIST1(J1),P(3*(J1-1)+1),P(3*(J1-1)+2),P(3*(J1-1)+3) 96:            WRITE(MYUNIT,'(2I6,3G20.10)') J1,NDIST1(J1),P(3*(J1-1)+1),P(3*(J1-1)+2),P(3*(J1-1)+3)
 97:         ENDDO 97:         ENDDO
 98:      ENDIF 98:      ENDIF
 99:   ENDIF 99:   ENDIF
100: 100: 
101:   DEALLOCATE(NDIST1,CON)101:   DEALLOCATE(NDIST1,CON)
102: 102: 
103: END SUBROUTINE PERC103: END SUBROUTINE PERC
104:  
105: SUBROUTINE PERCGROUP(P,NATOMS,PERCCUT,DEBUG,MYUNIT,RIGID) 
106: ! sf344> group the structures into connected parts, in order to move such parts together during MC step taking 
107: ! (should be useful for hierarchical global optimization where parts of the cluster should be moved  
108: ! around as a rigid body) 
109:  
110:   USE COMMONS, ONLY : atomingroup, modulart, modularcurrentn 
111:   IMPLICIT NONE 
112:  
113:   INTEGER, INTENT(IN)  :: NATOMS, MYUNIT 
114:   LOGICAL, INTENT(IN)  :: DEBUG, RIGID 
115:   DOUBLE PRECISION, INTENT(IN)  :: PERCCUT 
116:   INTEGER NSITES, J1, J2, NCYCLE, DMIN1, DMAX1, NUNCON1, groupindex, nsitesskip,nunconnsave,k1,k2 
117:   INTEGER, ALLOCATABLE :: NDIST1(:) 
118:   DOUBLE PRECISION DUMMY, COORDSNEW(3*NATOMS), P(3*NATOMS) !P is now intent inout 
119:   LOGICAL CHANGED, PERCT 
120:   LOGICAL, ALLOCATABLE :: CON(:,:) 
121:   
122:  
123:   IF(MODULART) THEN 
124:     NSITES = MODULARCURRENTN  
125:   ELSE 
126:     NSITES = NATOMS 
127:   END IF 
128:   IF (RIGID) THEN 
129:     NSITES = NSITES/2 
130:   ENDIF 
131:  
132: COORDSNEW(:)=P(:) 
133: P(:)=0.0D0  !reset P, we will populate it later in an ordered fashion 
134: NUNCONNSAVE=NSITES 
135: groupindex = 0 
136: ! WRITE(MYUNIT,'(A,2L5,3I6)') 'MODULART,RIGID,NSITES,NATOMS,MODULARCURRENTN=',MODULART,RIGID,NSITES,NATOMS,MODULARCURRENTN 
137:   IF(ALLOCATED(NDIST1)) DEALLOCATE(NDIST1) 
138:   IF(ALLOCATED(CON)) DEALLOCATE(CON) 
139:   ALLOCATE(NDIST1(NSITES), CON(NSITES,NSITES)) 
140:   ALLOCATE(atomingroup(nsites)) ! this will list the group index for each particle 
141: atomingroup(:)=0 
142: DO 
143: ! sf344>  after the first loop we have NUNNCONSAVE particles disconnected from particle 1, we need  
144: ! to check percolation only for these. Array P is being reordered so that each  
145: ! detected disconnected group is populating it sequentially. We take always AMERI...khm..CONNECTED particles first. 
146:  groupindex = groupindex + 1 
147:  NSITESSKIP=NSITES-NUNCONNSAVE !this is 0 for the first loop 
148:   IF(DEBUG) WRITE(MYUNIT,*) 'percgroup> Building group ', groupindex, nunconnsave, nsitesskip 
149: ! NSITES=NUNCONNSAVE 
150:  
151:   CON(1:NSITES,1:NSITES)=.FALSE. 
152:   DO J1=NSITESSKIP+1,NSITES 
153:     DO J2=J1+1,NSITES 
154:       DUMMY=(COORDSNEW(3*(J2-1)+1)-COORDSNEW(3*(J1-1)+1))**2+(COORDSNEW(3*(J2-1)+2)-COORDSNEW(3*(J1-1)+2))**2+(COORDSNEW(3*(J2-1)+3)-COORDSNEW(3*(J1-1)+3))**2 
155:       IF (DUMMY.LT.PERCCUT) THEN 
156:         CON(J2,J1)=.TRUE. 
157:         CON(J1,J2)=.TRUE. 
158:       ENDIF 
159:     ENDDO 
160:   ENDDO 
161:  
162: !  
163: ! Check that we have a percolating constraint network. 
164: ! 
165:   NDIST1(NSITESSKIP+1:NSITES)=1000000 
166:   NDIST1(NSITESSKIP+1)=0 
167:   NCYCLE=0 
168: 55 CHANGED=.FALSE. 
169:   NCYCLE=NCYCLE+1 
170:   DMIN1=100000 
171:   DMAX1=0 
172:   NUNCON1=0 
173:   DO J1=NSITESSKIP+1, NSITES 
174:     IF (NDIST1(J1).EQ.0) CYCLE ! first particle in the group 
175:     DO J2=NSITESSKIP+1, NSITES 
176:       IF (CON(J2,J1)) THEN 
177:         IF (NDIST1(J2)+1.LT.NDIST1(J1)) THEN 
178:           CHANGED=.TRUE. 
179:           NDIST1(J1)=NDIST1(J2)+1 
180:         ENDIF 
181:         IF (NDIST1(J1)+1.LT.NDIST1(J2)) THEN 
182:           CHANGED=.TRUE. 
183:           NDIST1(J2)=NDIST1(J1)+1 
184:         ENDIF 
185:       ENDIF 
186:     ENDDO 
187:     IF ((NDIST1(J1).GT.DMAX1).AND.(NDIST1(J1).NE.1000000)) DMAX1=NDIST1(J1) 
188:     IF (NDIST1(J1).LT.DMIN1) DMIN1=NDIST1(J1) 
189:     IF (NDIST1(J1).EQ.1000000) NUNCON1=NUNCON1+1 
190:   ENDDO 
191: !  PRINT *,'DMIN1,DMAX1,NUNCON1,NCYCLE,CHANGED=',DMIN1,DMAX1,NUNCON1,NCYCLE,CHANGED 
192:   IF (CHANGED) GOTO 55 
193:   IF (DEBUG) WRITE(MYUNIT,'(3(A,I8))') ' percgroup> steps to first atom in the group converged in ',NCYCLE-1, & 
194:      &                    ' cycles; maximum=',DMAX1,' disconnected=',NUNCON1 
195:   PERCT=.TRUE. 
196:   NUNCONNSAVE=NUNCON1 
197:   IF (NUNCON1.GT.0) THEN 
198: ! reorder the coordinates 
199:      K1=NSITESSKIP 
200:      K2=0 
201: !        WRITE(MYUNIT,*) 'before reordering', k1, k2, nsitesskip,nuncon1 
202: !        WRITE(MYUNIT,*) NDIST1(:) 
203: !        WRITE(MYUNIT,*) 'original coordinates: ' 
204: !        WRITE(MYUNIT,*) COORDSNEW(:) 
205:      DO J1=NSITESSKIP+1,NSITES 
206: !        WRITE(MYUNIT,*) 'sf344 debug> J1=', J1, NDIST1(J1) 
207:         IF(NDIST1(J1)<1000000) THEN 
208:           K1=K1+1 
209:           atomingroup(k1)=groupindex 
210:           P(3*K1-2:3*K1)=COORDSNEW(3*J1-2:3*J1) !copy over coordinates for connected group 
211:           IF(RIGID) P(3*K1-2+3*NATOMS/2:3*K1+3*NATOMS/2)=COORDSNEW(3*J1-2+3*NATOMS/2:3*J1+3*NATOMS/2) !copy over angle-axis coordinates as well 
212:         ELSE !disconnected particle, coordinates should go to the end of the array 
213:           K2=K2+1 
214: !          WRITE(MYUNIT,*) 'sf344 debug> ', 3*NSITES-3*NUNCON1+3*K2-2 
215: !          WRITE(MYUNIT,*) COORDSNEW(3*J1-2:3*J1) 
216:           P(3*NSITES-3*NUNCON1+3*K2-2:3*NSITES-3*NUNCON1+3*K2)=COORDSNEW(3*J1-2:3*J1) 
217:          IF(RIGID) P(3*NSITES-3*NUNCON1+3*K2+3*NATOMS/2-2:3*NSITES-3*NUNCON1+3*K2+3*NATOMS/2)=COORDSNEW(3*J1-2+3*NATOMS/2:3*J1+3*NATOMS/2) 
218:         END IF 
219:      END DO 
220: !        WRITE(MYUNIT,*) 'after reordering', k1, k2, nsitesskip 
221: !        WRITE(MYUNIT,*) 'reordered coordinates: ' 
222: !        WRITE(MYUNIT,*) P(:) 
223: ! save the reordered coordinate array 
224:      COORDSNEW(:)=P(:) 
225:      PERCT=.FALSE. 
226: !    IF (DEBUG) WRITE(MYUNIT,'(3G20.10)') P(1:3*NATOMS) 
227: !     IF (DEBUG) THEN !this debug printing does not make sense anymore, we already reordered the coordinates 
228: !        DO J1=1,NSITES 
229: !           WRITE(MYUNIT,'(2I6,3G20.10)') J1,NDIST1(J1),P(3*(J1-1)+1),P(3*(J1-1)+2),P(3*(J1-1)+3) 
230: !        ENDDO 
231: !     ENDIF 
232:   ELSE !NUNCON1==0, so every atom is connected in the last group, or the very last atom is disconnected 
233:    atomingroup(NSITESSKIP+1:NSITES)=groupindex 
234:      IF(DEBUG) THEN 
235:         WRITE(MYUNIT,*) 'percgroup> finished grouping, particles in the following groups: ' 
236:         WRITE(MYUNIT,*) atomingroup(:) 
237:      END IF 
238:      EXIT !exit the outer loop 
239:   ENDIF 
240: END DO 
241:  
242: !        WRITE(MYUNIT,*) 'percgroup> final reordered coordinates: ' 
243: !        WRITE(MYUNIT,*) P(:) 
244:  
245:   DEALLOCATE(NDIST1,CON) 
246: END SUBROUTINE PERCGROUP 


r31974/rigidmd.f90 2017-02-26 20:30:17.822913707 +0000 r31973/rigidmd.f90 2017-02-26 20:30:19.590937122 +0000
 54: ! RLSTSQ = RLIST**2 54: ! RLSTSQ = RLIST**2
 55: ! UPDATE = if true then update neighbor list 55: ! UPDATE = if true then update neighbor list
 56: ! PHI = an array keeping track of the angular deviation of each molecule 56: ! PHI = an array keeping track of the angular deviation of each molecule
 57: !       from it's initial position (at time t=NEQ).  It is calulated as 57: !       from it's initial position (at time t=NEQ).  It is calulated as
 58: !       a sum over time of W*DELT, so it is slightly approximate. 58: !       a sum over time of W*DELT, so it is slightly approximate.
 59: ! 59: !
 60: ! S = S is the additional degree of freedom in the Nose-Poincare thermostat 60: ! S = S is the additional degree of freedom in the Nose-Poincare thermostat
 61: ! PS = the conjugate momentum associated with S 61: ! PS = the conjugate momentum associated with S
 62: ! HNOT = the value of the hamiltonitan at time t=0.  Could also have been called 62: ! HNOT = the value of the hamiltonitan at time t=0.  Could also have been called
 63: !        HNAUGHT or H0   63: !        HNAUGHT or H0  
 64:       INTEGER, PARAMETER :: N = 1000, NN = N * 250 64:       INTEGER, PARAMETER :: N = 100, NN = N * 250
 65:       INTEGER          :: NMOL, NSITE=2, NTST, LIST(NN), POINT(N), G 65:       INTEGER          :: NMOL, NSITE=2, NTST, LIST(NN), POINT(N), G
 66:       DOUBLE PRECISION :: R(N,3), QTRN(N,4), V(N,3), W(N,3), F(N,3), T(N,3), T4(N,4), P(N,4), PHI(N,3)  66:       DOUBLE PRECISION :: R(N,3), QTRN(N,4), V(N,3), W(N,3), F(N,3), T(N,3), T4(N,4), P(N,4), PHI(N,3) 
 67:       DOUBLE PRECISION :: MST(3,3), R0(N,3)  67:       DOUBLE PRECISION :: MST(3,3), R0(N,3) 
 68:       DOUBLE PRECISION :: KEP(N), PEP(N), SUMEPT(N) 68:       DOUBLE PRECISION :: KEP(N), PEP(N), SUMEPT(N)
 69:       DOUBLE PRECISION :: BOXL, RCUT, RCUTSQ, RCUT2, RCUT2SQ, EPS4, CNSTA, CNSTB, DELT, HALFDT, MASS, JMI(3) 69:       DOUBLE PRECISION :: BOXL, RCUT, RCUTSQ, RCUT2, RCUT2SQ, EPS4, CNSTA, CNSTB, DELT, HALFDT, MASS, JMI(3)
 70:       DOUBLE PRECISION :: RLIST, RLSTSQ, TMPFIX 70:       DOUBLE PRECISION :: RLIST, RLSTSQ, TMPFIX
 71:       LOGICAL          :: UPDATE 71:       LOGICAL          :: UPDATE
 72:       !these are used only in the NVT RUN 72:       !these are used only in the NVT RUN
 73:       DOUBLE PRECISION :: FCTMQT=1, S=1, PS=1, HNOT=1 73:       DOUBLE PRECISION :: FCTMQT=1, S=1, PS=1, HNOT=1
 74:  74: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0