hdiff output

r31984/commons.f90 2017-02-28 10:30:07.025064477 +0000 r31983/commons.f90 2017-02-28 10:30:08.525084425 +0000
 66:      &                 CAPEPS2, CAPRAD, CAPRHO, CAPHEIGHT1, CAPHEIGHT2, & 66:      &                 CAPEPS2, CAPRAD, CAPRHO, CAPHEIGHT1, CAPHEIGHT2, &
 67:      &                 EPSR, GBKAPPA, GBKAPPRM, GBMU,GBNU, GBSIGNOT, GBEPSNOT, GBCHI, GBCHIPRM, & 67:      &                 EPSR, GBKAPPA, GBKAPPRM, GBMU,GBNU, GBSIGNOT, GBEPSNOT, GBCHI, GBCHIPRM, &
 68:      &                 SIGNOT, EPSNOT, SIGMAF, INVKAP, ESA(3), LPRSQ, LSQDFR, GBDPMU, GBDPEPS, GBDPFCT, & 68:      &                 SIGNOT, EPSNOT, SIGMAF, INVKAP, ESA(3), LPRSQ, LSQDFR, GBDPMU, GBDPEPS, GBDPFCT, &
 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), PERCGROUPSF, & 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, PERCGROUPCUT, &
 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, &
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, PERCGROUPT, &
110:      &        GENALT, MINDENSITYT, RESTRICTREGION, RESTRICTREGIONTEST, RESTRICTCYL, ACK1, ACK2, HARMONICF, PERCGROUPRESEEDT, &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, &
120:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &120:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &
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:       INTEGER, ALLOCATABLE :: atomingroup(:)
228:       DOUBLE PRECISION, ALLOCATABLE :: RESEEDCOORDS(:) 
229:       CHARACTER(LEN=80) :: MODULARRESTARTFILE 228:       CHARACTER(LEN=80) :: MODULARRESTARTFILE 
230: ! sf344> end of modular global optimization variables229: ! sf344> end of modular global optimization variables
231: 230: 
232:       INTEGER  NGBSITE, NPYSITE, PYBINARYTYPE1, PYSWAP(3), MAXINTERACTIONS231:       INTEGER  NGBSITE, NPYSITE, PYBINARYTYPE1, PYSWAP(3), MAXINTERACTIONS
233: 232: 
234:       DOUBLE PRECISION :: DZP1, DZP2, DZP3, DZP4, DZP5, DZP6, DZP7233:       DOUBLE PRECISION :: DZP1, DZP2, DZP3, DZP4, DZP5, DZP6, DZP7
235:       LOGICAL :: FIELDT, OHT, IHT, TDT, D5HT234:       LOGICAL :: FIELDT, OHT, IHT, TDT, D5HT
236:       DOUBLE PRECISION :: FOH, FIH, FTD, FD5H235:       DOUBLE PRECISION :: FOH, FIH, FTD, FD5H
237: 236: 
238:       !ds656> Box centroid to facilitate global237:       !ds656> Box centroid to facilitate global


r31984/gay-berne.f90 2017-02-28 10:30:07.277067828 +0000 r31983/gay-berne.f90 2017-02-28 10:30:08.777087780 +0000
4084:          E(3,1)  = -E(1,3)4084:          E(3,1)  = -E(1,3)
4085:          E(3,2)  = -E(2,3)4085:          E(3,2)  = -E(2,3)
4086:          E       = E/THETA4086:          E       = E/THETA
4087: 4087: 
4088:          ROTMAT = I3 + (1.D0-COS(THETA))*MATMUL(E,E) + E*SIN(THETA)4088:          ROTMAT = I3 + (1.D0-COS(THETA))*MATMUL(E,E) + E*SIN(THETA)
4089: 4089: 
4090:       ENDIF4090:       ENDIF
4091: 4091: 
4092:       END SUBROUTINE ELLIPSOIDROTATION4092:       END SUBROUTINE ELLIPSOIDROTATION
4093: 4093: 
 4094: 
 4095: 
4094: !     ----------------------------------------------------------------------------------------------4096: !     ----------------------------------------------------------------------------------------------
4095: 4097: 
4096:       SUBROUTINE TAKESTEPELPSD (NP)4098:       SUBROUTINE TAKESTEPELPSD (NP)
4097: 4099: 
4098: !     THIS ROUTINE TAKES STEP FOR SINGLE-SITE ELLIPSOIDAL BODIES ENSURING NO OVERLAP4100: !     THIS ROUTINE TAKES STEP FOR SINGLE-SITE ELLIPSOIDAL BODIES ENSURING NO OVERLAP
4099: 4101: 
4100:       USE COMMONS 4102:       USE COMMONS 
4101:       use pymodule, only : angle,angle2,twopi4103:       use pymodule, only : angle,angle2,twopi
4102:       use modamber9, only : ligtransstep, ligrotscale4104:       use modamber9, only : ligtransstep, ligrotscale
4103:       use rotations4105:       use rotations
4575: !                angle2=sqrt(angle2)4577: !                angle2=sqrt(angle2)
4576: !                angle = angle2 - int(angle2/twopi)*twopi4578: !                angle = angle2 - int(angle2/twopi)*twopi
4577: !                COORDS(J2-2:J2,NP)=COORDS(J2-2:J2,NP)/angle2 * angle4579: !                COORDS(J2-2:J2,NP)=COORDS(J2-2:J2,NP)/angle2 * angle
4578: !            end if4580: !            end if
4579: !        END DO4581: !        END DO
4580: 4582: 
4581: !RETURN4583: !RETURN
4582: 4584: 
4583:       END SUBROUTINE TAKESTEPELPSD4585:       END SUBROUTINE TAKESTEPELPSD
4584: 4586: 
4585:  
4586: !     ---------------------------------------------------------------------------------------------- 
4587:  
4588:       SUBROUTINE TAKESTEPPERCGROUPRESEED (NP) 
4589:  
4590: !     THIS ROUTINE TAKES STEP FOR SINGLE-SITE ELLIPSOIDAL BODIES ENSURING NO OVERLAP 
4591:  
4592:       USE COMMONS  
4593:       use pymodule, only : angle,angle2,twopi 
4594:       use modamber9, only : ligtransstep, ligrotscale 
4595:       use rotations 
4596:  
4597:       IMPLICIT NONE 
4598:  
4599:       INTEGER          :: NP, JMAX, JMAX2, REALNATOMS, OFFSET 
4600:       INTEGER          :: J1, J2, J3, J4, J5, J6 
4601: !      INTEGER          :: PTINDX, I, J, K 
4602:       LOGICAL          :: OVRLPT 
4603:  
4604:       DOUBLE PRECISION :: PI, DUMMY, DUMMY2, SAVECOORDS(3*NATOMS), Y(3*NATOMS) 
4605:       DOUBLE PRECISION :: DIST(3*NATOMS/2), XMASS, YMASS, ZMASS, DMAX, VMAX, VMAX2 
4606:       DOUBLE PRECISION :: VMIN, CMMAX, CMDIST(NATOMS/2), LOCALSTEP 
4607:       DOUBLE PRECISION :: DPRAND, RANDOM, THETA, PHI 
4608:       DOUBLE PRECISION :: AE(3,3), BE(3,3)  
4609:       DOUBLE PRECISION :: RI(3), RJ(3), RIJ(3) 
4610:       DOUBLE PRECISION :: P1(3), P2(3), ABSRIJ, RCUT, XMIN, FMIN, ECFVAL, ligandcoords(3), ligandcoordsrotated(3) 
4611: ! sf344 additions 
4612:       LOGICAL          :: DISSOCIATED(NATOMS/2),DISSOC 
4613:       INTEGER          :: CLOSESTATOMINDEX, K1, K2, i,j, k 
4614:       DOUBLE PRECISION :: MINDISTANCE, x1, x2, y1, y2, z1, z2, RCUTARRAY(NATOMS/2), RVEC(NATOMS/2,NATOMS/2) 
4615:       DOUBLE PRECISION, ALLOCATABLE :: atomblockcentre(:,:), rotationmatrix(:,:,:) 
4616:       double precision :: st,ct,sph,cph,sps,cosps, randomx, randomy, randomz, randomphi, randomtheta, randompsi 
4617:  
4618:       PI         = 4.D0*ATAN(1.0D0) 
4619:  
4620:       IF (PYGPERIODICT.OR.PYBINARYT) THEN 
4621:         DO J1=1,NATOMS/2 
4622:          RCUTARRAY(J1) = 2.0D0 * MAXVAL(PYA1bin(J1,:))  
4623:         END DO 
4624:          RCUT = MAXVAL(RCUTARRAY) 
4625: !        WRITE(*,*) 'RCUT=', RCUT 
4626:       ELSE IF (LJCAPSIDT) THEN 
4627:         PYA1bin(:,:)=0.5D0 
4628:         DO J1=1,NATOMS/2 
4629:          RCUTARRAY(J1) = 2.0D0 
4630:         END DO 
4631:          RCUT = 2.0D0 
4632:       ENDIF 
4633:  
4634:       REALNATOMS = NATOMS/2 
4635:       OFFSET     = 3 * REALNATOMS 
4636:  
4637:         SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,NP) 
4638: ! do the group moves if PERCGROUPT=.TRUE. 
4639:    IF(PERCGROUPT) THEN 
4640:         nblocks=atomingroup(REALNATOMS) !the last number in the atomingroup array is the total number of groups as well 
4641:         WRITE(MYUNIT,'(A,I6,A)') 'takesteppercgroupreseed> reseeding block centres from coords.reseed, and doing random rotation for ', nblocks, ' blocks'  
4642:         if(allocated(rotationmatrix)) deallocate(rotationmatrix) 
4643:         if(allocated(atomblockcentre)) deallocate(atomblockcentre) 
4644:         allocate(atomblockcentre(nblocks,3),rotationmatrix(nblocks,3,3)) 
4645: !WRITE(MYUNIT,*) 'before rotation and translation:' 
4646: !WRITE(MYUNIT,*) y(:) 
4647: OVRLPT=.TRUE. 
4648: 94  DO WHILE (OVRLPT) 
4649:  
4650:      y(:)=COORDS(:,NP) 
4651: ! random rotation matrices for each group 
4652:      IF(ligrotscale>0.0D0) then 
4653:        DO i=1,nblocks   
4654:         randomphi=(DPRAND()-0.5)*twopi*ligrotscale 
4655:         randomtheta=(DPRAND()-0.5)*pi*ligrotscale 
4656:         randompsi=(DPRAND()-0.5)*twopi*ligrotscale 
4657:         st=sin(randomtheta) 
4658:         ct=cos(randomtheta) 
4659:         sph=sin(randomphi) 
4660:         cph=cos(randomphi) 
4661:         sps=sin(randompsi) 
4662:         cosps=cos(randompsi) 
4663:  
4664:  
4665:         rotationmatrix(i,1,1)=cosps*cph-ct*sph*sps 
4666:         rotationmatrix(i,2,1)=cosps*sph+ct*cph*sps 
4667:         rotationmatrix(i,3,1)=sps*st 
4668:         rotationmatrix(i,1,2)=-sps*cph-ct*sph*cosps 
4669:         rotationmatrix(i,2,2)=-sps*sph+ct*cph*cosps 
4670:         rotationmatrix(i,3,2)=cosps*st 
4671:         rotationmatrix(i,1,3)=st*sph 
4672:         rotationmatrix(i,2,3)=-st*cph 
4673:         rotationmatrix(i,3,3)=ct 
4674:        END DO 
4675:  
4676:           atomblockcentre(:,:)=0.0D0 
4677:           do i=1,nblocks 
4678:              j=0 
4679:              do k=1,REALNATOMS 
4680:                if(i==atomingroup(k)) then 
4681:                 j=j+1 ! counter of atoms in each block 
4682:                 atomblockcentre(i,1)=atomblockcentre(i,1)+y(3*k-2) 
4683:                 atomblockcentre(i,2)=atomblockcentre(i,2)+y(3*k-1) 
4684:                 atomblockcentre(i,3)=atomblockcentre(i,3)+y(3*k  ) 
4685:                end if 
4686:              end do 
4687: !                write(MYUNIT,*) 'rotation matrix:', rotationmatrix(i,:,:) 
4688:                 atomblockcentre(i,:)=atomblockcentre(i,:)/j 
4689:                 atomblockcentre(i,:)=atomblockcentre(i,:) 
4690:           end do 
4691:           do i=1,nblocks 
4692:              do k=1,REALNATOMS 
4693:                if(i==atomingroup(k)) then 
4694: !                write(*,*) 'i, k=', i, k 
4695: ! move each block to the origin 
4696:                 ligandcoords(1)=y(3*k-2)-atomblockcentre(i,1) 
4697:                 ligandcoords(2)=y(3*k-1)-atomblockcentre(i,2) 
4698:                 ligandcoords(3)=y(3*k  )-atomblockcentre(i,3) 
4699: ! rotate the block of atoms with random rotation matrix 
4700:                 ligandcoordsrotated=MATMUL(rotationmatrix(i,:,:),ligandcoords) 
4701: !                WRITE(MYUNIT,*) 'ligandcoordsrotated: ', i, k, ligandcoordsrotated(:) 
4702: !                WRITE(MYUNIT,*) 'reseedcoords:', reseedcoords(:) 
4703: ! translate back the block of atoms to the original centre of their coordinates 
4704:                 y(3*k-2)=ligandcoordsrotated(1)+reseedcoords(3*i-2)*PERCGROUPSF 
4705:                 y(3*k-1)=ligandcoordsrotated(2)+reseedcoords(3*i-1)*PERCGROUPSF 
4706:                 y(3*k  )=ligandcoordsrotated(3)+reseedcoords(3*i  )*PERCGROUPSF 
4707:                 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,:,:))) 
4708:                end if 
4709:              end do 
4710:           end do 
4711: !                WRITE(MYUNIT,*) 'after rotation:' 
4712: !                WRITE(MYUNIT,*) y(:) 
4713:      END IF !ligrotscale 
4714:  
4715:       COORDS(:,NP)=y(:) 
4716:  
4717: ! check for overlap 
4718:       DO J1 = 1, REALNATOMS 
4719:  
4720:          J3 = 3*J1 
4721:          J5 = OFFSET + J3 
4722:  
4723:             OVRLPT = .FALSE. 
4724:  
4725:             RI(:)  = COORDS(J3-2:J3,NP) 
4726:             P1(:)  = COORDS(J5-2:J5,NP) 
4727:  
4728: !     ROTATION MATRIX 
4729:  
4730:             CALL ELPSDRTMT (P1, PYA1bin(J1,:), AE) 
4731:  
4732:             DO J2 = 1, REALNATOMS 
4733:  
4734:                IF (J2 == J1) CYCLE 
4735:   
4736:                J4    = 3*J2 
4737:                J6    = OFFSET + J4 
4738:                RJ(:) = COORDS(J4-2:J4,NP) 
4739:                P2(:) = COORDS(J6-2:J6,NP) 
4740:  
4741: !     ROTATION MATRIX 
4742:  
4743:                CALL ELPSDRTMT (P2, PYA1bin(J2,:), BE) 
4744:  
4745:                RIJ    = RI - RJ 
4746:                ABSRIJ = DSQRT(DOT_PRODUCT(RIJ,RIJ)) 
4747:                          
4748:                IF (ABSRIJ < RCUT) THEN  
4749:  
4750: !     DETERMINE ELLIPTIC CONTACT FUNCTION 
4751:  
4752:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE, BE, RIJ, XMIN, FMIN) 
4753:  
4754:                   ECFVAL = - FMIN 
4755: ! allow for a slight overlap  
4756:                   IF (ECFVAL < PYOVERLAPTHRESH) THEN 
4757:                    OVRLPT = .TRUE. 
4758:                    PERCGROUPSF=PERCGROUPSF*1.1D0 
4759:                    IF(DEBUG) WRITE(MYUNIT,*) 'takestepelpsd> PERCGROUP - particles overlapping, scaling up the reseeded centres' 
4760:                    IF(DEBUG) WRITE(MYUNIT,*) 'takestepelpsd> scale factor, J1, J2, ECF', PERCGROUPSF,J1, J2, ECFVAL 
4761:                    GO TO 94 
4762:                   ENDIF 
4763:  
4764:                ENDIF 
4765:  
4766:             ENDDO  ! END LOOP OVER J2 
4767:       ENDDO  ! END LOOP OVER J1    
4768:  
4769:     END DO ! DO WHILE(OVERLAPT) 
4770:    END IF ! PERCGROUPT 
4771: !WRITE(MYUNIT,*) 'rotated and translated coordinates:' 
4772: !WRITE(MYUNIT,*) y(:) 
4773:  
4774:  
4775:       IF (FREEZE) THEN 
4776:          DO J2=1,NATOMS 
4777:            IF (FROZEN(J2)) THEN 
4778:               COORDS(3*(J2-1)+1:3*(J2-1)+3,NP)=SAVECOORDS(3*(J2-1)+1:3*(J2-1)+3) 
4779:            ENDIF 
4780:          ENDDO 
4781:       ENDIF 
4782:  
4783: RETURN 
4784:  
4785:       END SUBROUTINE TAKESTEPPERCGROUPRESEED 
4786:  
4787: SUBROUTINE TAKESTEPSWAPMOVES(NP)4587: SUBROUTINE TAKESTEPSWAPMOVES(NP)
4788: use commons, only : NATOMS, COORDS, PYBINARYTYPE1,MYUNIT,PYSWAP!,SWAPMOVEST4588: use commons, only : NATOMS, COORDS, PYBINARYTYPE1,MYUNIT,PYSWAP!,SWAPMOVEST
4789: 4589: 
4790: implicit none4590: implicit none
4791: 4591: 
4792: integer :: NP, RANDOM1, RANDOM2, J1, J24592: integer :: NP, RANDOM1, RANDOM2, J1, J2
4793: double precision :: COORDSSTORE(3,4),DPRAND4593: double precision :: COORDSSTORE(3,4),DPRAND
4794: 4594: 
4795: ! get two random integer numbers, one for atom index type 1, one for atom index type 24595: ! get two random integer numbers, one for atom index type 1, one for atom index type 2
4796: 4596: 


r31984/keywords.f 2017-02-28 10:30:07.529071178 +0000 r31983/keywords.f 2017-02-28 10:30:09.429096571 +0000
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.0D1001005:       PERCGROUPCUT=1.0D100
1006:       PERCGROUPT=.FALSE.1006:       PERCGROUPT=.FALSE.
1007:       PERCGROUPRESEEDT=.FALSE. 
1008:       PERCGROUPSF=1.0D0 
1009:       GENALT=.FALSE.1007:       GENALT=.FALSE.
1010: !1008: !
1011: ! Reservoir list parameters for PT / BSPT.1009: ! Reservoir list parameters for PT / BSPT.
1012: !1010: !
1013:       RESERVOIRT=.FALSE.1011:       RESERVOIRT=.FALSE.
1014:       USERES=01012:       USERES=0
1015:       EXEQ=01013:       EXEQ=0
1016:       PERTSTEP=0.01D01014:       PERTSTEP=0.01D0
1017:       RES_PSWAP = 0.0D01015:       RES_PSWAP = 0.0D0
1018:       BETA_RES = 0.0D01016:       BETA_RES = 0.0D0
5171:          PERCCOMPMARKOV=.FALSE.5169:          PERCCOMPMARKOV=.FALSE.
5172:          IF (NITEMS.GT.1) CALL READF(PERCCUT)5170:          IF (NITEMS.GT.1) CALL READF(PERCCUT)
5173:          IF (NITEMS.GT.2) CALL READF(K_COMP)5171:          IF (NITEMS.GT.2) CALL READF(K_COMP)
5174:          IF (NITEMS.GT.3) CALL READF(GUIDECUT)5172:          IF (NITEMS.GT.3) CALL READF(GUIDECUT)
5175: 5173: 
5176:       ELSE IF (WORD.EQ.'PERCGROUP') THEN5174:       ELSE IF (WORD.EQ.'PERCGROUP') THEN
5177:          PERCGROUPT=.TRUE.5175:          PERCGROUPT=.TRUE.
5178:          CALL READF(PERCGROUPCUT)5176:          CALL READF(PERCGROUPCUT)
5179:          CALL READF(LIGTRANSSTEP)5177:          CALL READF(LIGTRANSSTEP)
5180:          CALL READF(LIGROTSCALE)5178:          CALL READF(LIGROTSCALE)
5181:          IF(NITEMS.GT.4) THEN !seed groups with predefined coordinates 
5182:              PERCGROUPRESEEDT=.TRUE. 
5183:              CALL READF(PERCGROUPSF) 
5184:          END IF 
5185: 5179: 
5186:       ELSE IF (WORD.EQ.'PERIODIC') THEN5180:       ELSE IF (WORD.EQ.'PERIODIC') THEN
5187:          PERIODIC=.TRUE.5181:          PERIODIC=.TRUE.
5188:          CALL READF(XX)5182:          CALL READF(XX)
5189:          BOXLX=XX5183:          BOXLX=XX
5190:          BOX3D(1) = XX5184:          BOX3D(1) = XX
5191:          IF (NITEMS.GT.2) THEN5185:          IF (NITEMS.GT.2) THEN
5192:             CALL READF(XX)5186:             CALL READF(XX)
5193:             BOXLY=XX5187:             BOXLY=XX
5194:             BOX3D(2) = XX5188:             BOX3D(2) = XX


r31984/mc.F 2017-02-28 10:30:07.781074535 +0000 r31983/mc.F 2017-02-28 10:30:09.685099857 +0000
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) THEN1002:                 IF(PERCGROUPT) THEN
1003: ! sf344> group particles according to percolation, reorder them1003: ! sf344> group particles according to percolation, reorder them
1004:                   CALL PERCGROUP(COORDS(1:3*NATOMS,JP),NATOMS,PERCGROUPCUT,DEBUG,MYUNIT,RIGID)1004:                   CALL PERCGROUP(COORDS(1:3*NATOMS,JP),NATOMS,PERCGROUPCUT,DEBUG,MYUNIT,RIGID)
1005:                   IF(NQ(JP)==0.AND.PERCGROUPRESEEDT) THEN !reseeding with coordinates in file coords.reseed1005: !                 in this case COORDSO is the same as COORDS, except atoms are being reshuffled. Reset
1006:                     OPEN(unit=299,file='coords.reseed',status='old')1006: !                  COORDSO to get rid of sanity check later on
1007:                     IF(ALLOCATED(RESEEDCOORDS)) DEALLOCATE(RESEEDCOORDS)1007:                   COORDSO(:,JP)=COORDS(:,JP)
1008:                     ALLOCATE(RESEEDCOORDS(3*atomingroup(NATOMS/2))) 
1009:                     DO j1=1,atomingroup(NATOMS/2) ! this is the total number of groups 
1010:                         ! we have to have just as many coordinates in coords.reseed 
1011:                        READ(299,*) RESEEDCOORDS(3*j1-2), RESEEDCOORDS(3*j1-1), RESEEDCOORDS(3*j1)   
1012:                     END DO 
1013:                     CLOSE(299) 
1014:                     CALL TAKESTEPPERCGROUPRESEED(JP) 
1015:                   END IF 
1016:                 END IF1008:                 END IF
1017: ! sf344> particles are now ordered according to the percolated groups they are in. Group indices are1009: ! sf344> particles are now ordered according to the percolated groups they are in. Group indices are
1018: !        saved in atomingroup(:)1010: !        saved in atomingroup(:)
1019: !               seed the random number generator with system time + MYNODE (for MPI runs)1011: !               seed the random number generator with system time + MYNODE (for MPI runs)
1020:                   IF(RANDOMSEEDT) THEN1012:                   IF(RANDOMSEEDT) THEN
1021:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)1013:                      CALL DATE_AND_TIME(datechar,timechar,zonechar,values)
1022:                      itime1= values(6)*60 + values(7)1014:                      itime1= values(6)*60 + values(7)
1023:                      CALL SDPRND(itime1+MYNODE)1015:                      CALL SDPRND(itime1+MYNODE)
1024:                   END IF1016:                   END IF
1025:                   IF(SWAPMOVEST) THEN1017:                   IF(SWAPMOVEST) THEN
1026:                      CALL TAKESTEPSWAPMOVES(JP)1018:                      CALL TAKESTEPSWAPMOVES(JP)
1027:                      CALL TAKESTEPELPSD(JP)1019:                      CALL TAKESTEPELPSD(JP)
1028:                   ELSE1020:                   ELSE
1029:                    ! avoid taking steps for the first quench if we are reseeding the group centres1021:                      CALL TAKESTEPELPSD(JP)
1030:                     IF((NQ(JP)>0.AND.PERCGROUPRESEEDT).OR.(.NOT.PERCGROUPRESEEDT)) THEN  
1031:                         CALL TAKESTEPELPSD(JP) 
1032:                     END IF 
1033:                   END IF1022:                   END IF
1034: 1023: 
1035: !1024: !
1036: ! MARK'S MIXED CLUSTER STEP TAKING1025: ! MARK'S MIXED CLUSTER STEP TAKING
1037: !1026: !
1038: !     MAM1000 >  Charged/neutral particle swaps for LJ+Coulomb mixed clusters1027: !     MAM1000 >  Charged/neutral particle swaps for LJ+Coulomb mixed clusters
1039:                ELSE IF (LJCOULT) THEN1028:                ELSE IF (LJCOULT) THEN
1040:                   RANDOM=DPRAND()1029:                   RANDOM=DPRAND()
1041:                   IF (RANDOM < COULSWAP) THEN1030:                   IF (RANDOM < COULSWAP) THEN
1042:                      CALL TAKESTEPLJC(JP)1031:                      CALL TAKESTEPLJC(JP)
1967:                         EPREV(JP)=HBONDESAVE1956:                         EPREV(JP)=HBONDESAVE
1968: !                        WRITE(MYUNIT,'(A)') " HBONDMATRIX> Markov energy reset to pass sanity checks"1957: !                        WRITE(MYUNIT,'(A)') " HBONDMATRIX> Markov energy reset to pass sanity checks"
1969:                      ENDIF1958:                      ENDIF
1970:                   ENDIF1959:                   ENDIF
1971:                ENDIF1960:                ENDIF
1972: !1961: !
1973: !  Sanity check to make sure the Markov energy agrees with COORDSO. 1962: !  Sanity check to make sure the Markov energy agrees with COORDSO. 
1974: !  Stop if not true.1963: !  Stop if not true.
1975: !1964: !
1976:                IF ((DEBUG.OR.CHECKMARKOVT).AND.(.NOT.(HBONDMATRIX.OR.PERMOPT.OR.PERMINVOPT.OR.FEBHT.OR.QALCST.OR. 1965:                IF ((DEBUG.OR.CHECKMARKOVT).AND.(.NOT.(HBONDMATRIX.OR.PERMOPT.OR.PERMINVOPT.OR.FEBHT.OR.QALCST.OR. 
1977:      &                           MODULART.OR.PERCGROUPT))) THEN1966: &                               MODULART.OR.PERCGROUPT))) THEN
1978: ! jwrm2> If we're using percolate, we need to set compression to whether it was on or not when calculating the Markov energy1967: ! jwrm2> If we're using percolate, we need to set compression to whether it was on or not when calculating the Markov energy
1979:                   QUENCHCOMP = COMPON1968:                   QUENCHCOMP = COMPON
1980:                   IF (PERCOLATET) COMPON = PERCCOMPMARKOV1969:                   IF (PERCOLATET) COMPON = PERCCOMPMARKOV
1981:                   CALL POTENTIAL(COORDSO(:,JP),GRAD,OPOTEL,.FALSE.,.FALSE.)1970:                   CALL POTENTIAL(COORDSO(:,JP),GRAD,OPOTEL,.FALSE.,.FALSE.)
1982:                   IF (ABS(OPOTEL-EPREV(JP)).GT.ABS(ECONV)) THEN1971:                   IF (ABS(OPOTEL-EPREV(JP)).GT.ABS(ECONV)) THEN
1983:                      IF (EVAP) THEN1972:                      IF (EVAP) THEN
1984:                         WRITE(MYUNIT,'(3(A,G20.10))') 'mc> WARNING - energy for saved coordinates ',OPOTEL,1973:                         WRITE(MYUNIT,'(3(A,G20.10))') 'mc> WARNING - energy for saved coordinates ',OPOTEL,
1985:      &                     ' differs from Markov energy ',EPREV(JP),' because an atom moved outside the container'1974:      &                     ' differs from Markov energy ',EPREV(JP),' because an atom moved outside the container'
1986:                      ELSE1975:                      ELSE
1987:                         WRITE(MYUNIT,'(2(A,G20.10))') 'mc> ERROR - energy for coordinates in COORDSO=',OPOTEL,1976:                         WRITE(MYUNIT,'(2(A,G20.10))') 'mc> ERROR - energy for coordinates in COORDSO=',OPOTEL,


r31984/perc.f90 2017-02-28 10:30:08.029077834 +0000 r31983/perc.f90 2017-02-28 10:30:09.933103155 +0000
130:   ENDIF130:   ENDIF
131: 131: 
132: COORDSNEW(:)=P(:)132: COORDSNEW(:)=P(:)
133: P(:)=0.0D0  !reset P, we will populate it later in an ordered fashion133: P(:)=0.0D0  !reset P, we will populate it later in an ordered fashion
134: NUNCONNSAVE=NSITES134: NUNCONNSAVE=NSITES
135: groupindex = 0135: groupindex = 0
136: ! WRITE(MYUNIT,'(A,2L5,3I6)') 'MODULART,RIGID,NSITES,NATOMS,MODULARCURRENTN=',MODULART,RIGID,NSITES,NATOMS,MODULARCURRENTN136: ! WRITE(MYUNIT,'(A,2L5,3I6)') 'MODULART,RIGID,NSITES,NATOMS,MODULARCURRENTN=',MODULART,RIGID,NSITES,NATOMS,MODULARCURRENTN
137:   IF(ALLOCATED(NDIST1)) DEALLOCATE(NDIST1)137:   IF(ALLOCATED(NDIST1)) DEALLOCATE(NDIST1)
138:   IF(ALLOCATED(CON)) DEALLOCATE(CON)138:   IF(ALLOCATED(CON)) DEALLOCATE(CON)
139:   ALLOCATE(NDIST1(NSITES), CON(NSITES,NSITES))139:   ALLOCATE(NDIST1(NSITES), CON(NSITES,NSITES))
140:   IF(ALLOCATED(atomingroup)) DEALLOCATE(atomingroup) 
141:   ALLOCATE(atomingroup(nsites)) ! this will list the group index for each particle140:   ALLOCATE(atomingroup(nsites)) ! this will list the group index for each particle
142: atomingroup(:)=0141: atomingroup(:)=0
143: DO142: DO
144: ! sf344>  after the first loop we have NUNNCONSAVE particles disconnected from particle 1, we need 143: ! sf344>  after the first loop we have NUNNCONSAVE particles disconnected from particle 1, we need 
145: ! to check percolation only for these. Array P is being reordered so that each 144: ! to check percolation only for these. Array P is being reordered so that each 
146: ! detected disconnected group is populating it sequentially. We take always AMERI...khm..CONNECTED particles first.145: ! detected disconnected group is populating it sequentially. We take always AMERI...khm..CONNECTED particles first.
147:  groupindex = groupindex + 1146:  groupindex = groupindex + 1
148:  NSITESSKIP=NSITES-NUNCONNSAVE !this is 0 for the first loop147:  NSITESSKIP=NSITES-NUNCONNSAVE !this is 0 for the first loop
149:   IF(DEBUG) WRITE(MYUNIT,*) 'percgroup> Building group ', groupindex, nunconnsave, nsitesskip148:   IF(DEBUG) WRITE(MYUNIT,*) 'percgroup> Building group ', groupindex, nunconnsave, nsitesskip
150: ! NSITES=NUNCONNSAVE149: ! NSITES=NUNCONNSAVE


r31984/sandbox.f90 2017-02-28 10:30:08.277081133 +0000 r31983/sandbox.f90 2017-02-28 10:30:10.181106451 +0000
779:     end do779:     end do
780: 780: 
781:     ! initialize molecules781:     ! initialize molecules
782:     do mol_index = 1, size(molecules)782:     do mol_index = 1, size(molecules)
783:         call initialize_sandbox_molecule(molecules(mol_index), mol_index)783:         call initialize_sandbox_molecule(molecules(mol_index), mol_index)
784:     end do784:     end do
785: 785: 
786: end subroutine sandbox_input786: end subroutine sandbox_input
787: 787: 
788: subroutine sandbox(x, grad, energy, gtest)788: subroutine sandbox(x, grad, energy, gtest)
789:     use commons, only: vt, twod, frozen789:     use commons, only: vt, twod
790:     use sandbox_module790:     use sandbox_module
791:     implicit none791:     implicit none
792:     double precision, intent(inout) :: x(x_size)792:     double precision, intent(inout) :: x(x_size)
793:     logical, intent(in) :: gtest793:     logical, intent(in) :: gtest
794:     double precision, intent(out) :: energy, grad(x_size)794:     double precision, intent(out) :: energy, grad(x_size)
795: 795: 
796:     type(sandbox_molecule), pointer :: moli, molj796:     type(sandbox_molecule), pointer :: moli, molj
797:     type(sandbox_site), pointer :: sitei, sitej797:     type(sandbox_site), pointer :: sitei, sitej
798:     integer moli_index, molj_index, sitei_index, sitej_index798:     integer moli_index, molj_index, sitei_index, sitej_index
799:     double precision :: energy_contrib, grad_contrib(12)799:     double precision :: energy_contrib, grad_contrib(12)
847:             end do847:             end do
848:         end do848:         end do
849:     end do849:     end do
850: 850: 
851:     ! 2D: set all z gradients to zero851:     ! 2D: set all z gradients to zero
852:     if(twod) then852:     if(twod) then
853:         do moli_index = 1, num_mols853:         do moli_index = 1, num_mols
854:             grad(molecules(moli_index)%r_index + 2) = 0.D0854:             grad(molecules(moli_index)%r_index + 2) = 0.D0
855:         end do855:         end do
856:     end if856:     end if
857:     ! freeze keyword: set gradients to zero 
858:     if(gtest .and. any(frozen)) then 
859: !        STOP 
860:         do moli_index = 1, num_mols 
861:             moli => molecules(moli_index) 
862:             if(frozen(moli_index)) grad(moli%r_index: moli%r_index + 2) = 0.d0 
863:             if(frozen(moli_index + num_mols)) grad(moli%p_index: moli%p_index + 2) = 0.d0 
864:         end do 
865:     end if 
866: 857: 
867: end subroutine sandbox858: end subroutine sandbox
868: 859: 
869: subroutine sandbox_output860: subroutine sandbox_output
870:     use commons, only: natoms, nsave, mpit, mynode861:     use commons, only: natoms, nsave, mpit, mynode
871:     use qmodule, only: qmin, qminp, ff862:     use qmodule, only: qmin, qminp, ff
872:     use sandbox_module863:     use sandbox_module
873:     implicit none864:     implicit none
874: 865: 
875:     integer :: sandout_unit, coords_unit866:     integer :: sandout_unit, coords_unit
931:         end do922:         end do
932:         close(coords_unit)923:         close(coords_unit)
933: 924: 
934:     end do925:     end do
935: 926: 
936:     close(sandout_unit)927:     close(sandout_unit)
937:             928:             
938: end subroutine sandbox_output929: end subroutine sandbox_output
939: 930: 
940: subroutine sandbox_takestep(np)931: subroutine sandbox_takestep(np)
941:     use commons, only: myunit, coords, radius, percolatet, perccut, tmove, omove, step, ostep, astep, vt, twod, frozen932:     use commons, only: myunit, coords, radius, percolatet, perccut, tmove, omove, step, ostep, astep, vt, twod
942:     use sandbox_module933:     use sandbox_module
943:     use vec3, only: vec_len, vec_random934:     use vec3, only: vec_len, vec_random
944:     use rotations, only: rot_takestep_aa935:     use rotations, only: rot_takestep_aa
945: 936: 
946:     implicit none937:     implicit none
947:     integer, intent(in) :: np938:     integer, intent(in) :: np
948: 939: 
949:     double precision, target :: x(size(coords(:, np)))940:     double precision, target :: x(size(coords(:, np)))
950:     double precision, pointer :: r(:), p(:)941:     double precision, pointer :: r(:), p(:)
951:     double precision :: step_size, min_vt, dice_roll, twod_step(3)942:     double precision :: step_size, min_vt, dice_roll, twod_step(3)
980:         if(twod) then971:         if(twod) then
981:             ! translational step972:             ! translational step
982:             do973:             do
983:                 call random_number(twod_step(1))974:                 call random_number(twod_step(1))
984:                 call random_number(twod_step(2))975:                 call random_number(twod_step(2))
985:                 twod_step(3) = 0.D0976:                 twod_step(3) = 0.D0
986:                 if(vec_len(twod_step) <= 1.0) exit977:                 if(vec_len(twod_step) <= 1.0) exit
987:             end do978:             end do
988:             r(:) = r(:) + (step(np) * twod_step(:))979:             r(:) = r(:) + (step(np) * twod_step(:))
989:         else980:         else
990:          if(.not.(frozen(moli_index) .or. frozen(moli_index + num_mols))) then 
991:             ! angular move: accept step with probability 1 - exp(-beta * delta_energy), where beta = 1/kT = astep981:             ! angular move: accept step with probability 1 - exp(-beta * delta_energy), where beta = 1/kT = astep
992:             if(astep(np) > 0.D0) then982:             if(astep(np) > 0.D0) then
993:                 call random_number(dice_roll)983:                 call random_number(dice_roll)
994:                 if(dice_roll > exp(-astep(np) * (vt(moli_index) - min_vt))) then984:                 if(dice_roll > exp(-astep(np) * (vt(moli_index) - min_vt))) then
995:                     r = max_distance(molecules) * vec_random()985:                     r = max_distance(molecules) * vec_random()
996: 986: 
997:                     ! if using percolate, then pull the molecule in until it is within percolation distance of another987:                     ! if using percolate, then pull the molecule in until it is within percolation distance of another
998:                     if(percolatet) then988:                     if(percolatet) then
999:                         stop_pulling = .false.989:                         stop_pulling = .false.
1000:                         do990:                         do
1013:                         end do1003:                         end do
1014:                     else ! if using radius, pull in until the molecule is inside the container1004:                     else ! if using radius, pull in until the molecule is inside the container
1015:                         do1005:                         do
1016:                             if(vec_len(r) < radius) exit1006:                             if(vec_len(r) < radius) exit
1017:                             r(:) = 0.95 * r(:)1007:                             r(:) = 0.95 * r(:)
1018:                         end do1008:                         end do
1019:                     end if1009:                     end if
1020: 1010: 
1021:                 end if1011:                 end if
1022:             end if1012:             end if
1023:          end if 1013: 
1024:             ! translational move: uniformly distributed in a sphere of radius step(np)1014:             ! translational move: uniformly distributed in a sphere of radius step(np)
1025:          if(.not. frozen(moli_index)) then 
1026:             if(tmove(np)) then1015:             if(tmove(np)) then
1027:                 call random_number(step_size)1016:                 call random_number(step_size)
1028:                 step_size = sqrt(step(np) * step_size)  ! need to take square root to get uniform volume distribution1017:                 step_size = sqrt(step(np) * step_size)  ! need to take square root to get uniform volume distribution
1029:                 r(:) = r(:) + step_size * vec_random()1018:                 r(:) = r(:) + step_size * vec_random()
1030:             end if1019:             end if
1031:          end if 
1032:         end if ! 2d1020:         end if ! 2d
1033: 1021: 
1034:         ! orientational move1022:         ! orientational move
1035:          if(.not. frozen(moli_index + num_mols)) then1023:         if(omove(np)) then
1036:             if(omove(np)) then1024:             call random_number(step_size)
1037:               call random_number(step_size)1025:             call rot_takestep_aa(p(:), step_size * ostep(np))
1038:               call rot_takestep_aa(p(:), step_size * ostep(np))1026:         end if
1039:             end if1027: 
1040:          end if 
1041:     end do1028:     end do
1042: 1029: 
1043:     ! copy local results back1030:     ! copy local results back
1044:     coords(:, np) = x(:)1031:     coords(:, np) = x(:)
1045: 1032: 
1046: end subroutine sandbox_takestep1033: end subroutine sandbox_takestep


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0