hdiff output

r30520/geopt.f 2016-05-26 15:30:14.088123374 +0100 r30519/geopt.f 2016-05-26 15:30:16.316154455 +0100
390:             ENDIF390:             ENDIF
391:             IF (RINGPOLYMERT) THEN391:             IF (RINGPOLYMERT) THEN
392:                IF (RPSYSTEM(1:4).EQ.'AECK') THEN392:                IF (RPSYSTEM(1:4).EQ.'AECK') THEN
393:                   NEXMODES=0393:                   NEXMODES=0
394:                ELSE394:                ELSE
395:                   NEXMODES=6395:                   NEXMODES=6
396:                ENDIF396:                ENDIF
397:                IF (BFGSTST.OR.(INR.EQ.2)) NEXMODES=NEXMODES+1397:                IF (BFGSTST.OR.(INR.EQ.2)) NEXMODES=NEXMODES+1
398:                IF (BFGSMINT) NEXMODES=NEXMODES+2398:                IF (BFGSMINT) NEXMODES=NEXMODES+2
399:             ENDIF399:             ENDIF
400:             IF (STEALTHYT) THEN 
401:                IF (ENERGY.LT.1.0D-8) THEN 
402:                   NEXMODES=NATOMS*3 - STM*2 
403:                ELSE 
404:                   NEXMODES=3 
405:                ENDIF 
406:             ENDIF 
407:             IF (VARIABLES.OR.MIEFT) NEXMODES=0400:             IF (VARIABLES.OR.MIEFT) NEXMODES=0
408:             WRITE(*,'(A,I6)') ' geopt> Number of zero/imaginary eigenvalues assumed to be ',NEXMODES401:             WRITE(*,'(A,I6)') ' geopt> Number of zero/imaginary eigenvalues assumed to be ',NEXMODES
409:  
410:             IF (LOWESTFRQT) THEN402:             IF (LOWESTFRQT) THEN
411: C403: C
412: C  Calculate lowest non-zero eigenvalue and dump to min.data.info file404: C  Calculate lowest non-zero eigenvalue and dump to min.data.info file
413: C  No mass-weighting here!405: C  No mass-weighting here!
414: C  'U' specifies that the upper triangle contains the Hessian.406: C  'U' specifies that the upper triangle contains the Hessian.
415: C  Need to save HESS before this call and restore afterwards.407: C  Need to save HESS before this call and restore afterwards.
416: C408: C
417:     409:     
418:                ABSTOL=DLAMCH('Safe  minimum')410:                ABSTOL=DLAMCH('Safe  minimum')
419:                IF (ALLOCATED(ZWK)) DEALLOCATE(ZWK)411:                IF (ALLOCATED(ZWK)) DEALLOCATE(ZWK)
877:                            RBAANORMALMODET = .FALSE.869:                            RBAANORMALMODET = .FALSE.
878:                         ELSE870:                         ELSE
879:                            CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NENDHESS,ABSTOL,871:                            CALL DSYEVR('N','I','U',NOPT,HESS,NOPT,0.0D0,1.0D0,1,NENDHESS,ABSTOL,
880:      &                          NFOUND,EVALUES,ZWK,NOPT,ISUPPZ,WORK,LWORK,IWORK,ILWORK,INFO )872:      &                          NFOUND,EVALUES,ZWK,NOPT,ISUPPZ,WORK,LWORK,IWORK,ILWORK,INFO )
881:                         ENDIF873:                         ENDIF
882:                         DEALLOCATE(ZWK)874:                         DEALLOCATE(ZWK)
883:                      ENDIF875:                      ENDIF
884:                   ENDIF876:                   ENDIF
885:                ENDIF877:                ENDIF
886:             ENDIF878:             ENDIF
887:  
888: C879: C
889: C  MASSWT2 and MASSWTRP do not mass weight Q and VNEW, but MASSWT does. Need to undo this880: C  MASSWT2 and MASSWTRP do not mass weight Q and VNEW, but MASSWT does. Need to undo this
890: C  if DUMPDATAT is .TRUE. for comparison with pathsample (which uses unit masses).881: C  if DUMPDATAT is .TRUE. for comparison with pathsample (which uses unit masses).
891: C  Probably best to always undo it, in case we need non-mass-weighted Q somewhere882: C  Probably best to always undo it, in case we need non-mass-weighted Q somewhere
892: C  further down.883: C  further down.
893: C884: C
894: C           IF (DUMPDATAT.AND.(.NOT.(CHRMMT.OR.RINGPOLYMERT))) THEN885: C           IF (DUMPDATAT.AND.(.NOT.(CHRMMT.OR.RINGPOLYMERT))) THEN
895:             IF (.NOT. (CHRMMT .OR. RBAAT)) THEN ! hk286886:             IF (.NOT. (CHRMMT .OR. RBAAT)) THEN ! hk286
896:                IF (.NOT.(RIGIDINIT.OR.VARIABLES)) THEN ! jmc49 bugfix, because the mass weighting is not done if rigidinit 887:                IF (.NOT.(RIGIDINIT.OR.VARIABLES)) THEN ! jmc49 bugfix, because the mass weighting is not done if rigidinit 
897:                   DO J1=1,NATOMS888:                   DO J1=1,NATOMS


r30520/grad.f90 2016-05-26 15:30:13.412113943 +0100 r30519/grad.f90 2016-05-26 15:30:15.652145192 +0100
239:   &               PARAM1*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1))/PARAM1))**2+ &239:   &               PARAM1*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1))/PARAM1))**2+ &
240:   &                           (XYZ(NOPT*(J1+1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2) - &240:   &                           (XYZ(NOPT*(J1+1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2) - &
241:   &               PARAM2*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2))/PARAM2))**2241:   &               PARAM2*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2))/PARAM2))**2
242:    IF (.NOT.TWOD) DISTA=DISTA+(XYZ(NOPT*(J1+1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3) - &242:    IF (.NOT.TWOD) DISTA=DISTA+(XYZ(NOPT*(J1+1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3) - &
243:   &               PARAM3*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3))/PARAM3))**2243:   &               PARAM3*NINT((XYZ(NOPT*(J1+1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3))/PARAM3))**2
244:                         DISTB=DISTB+ &244:                         DISTB=DISTB+ &
245:   &                           (XYZ(NOPT*(J1-1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1) - &245:   &                           (XYZ(NOPT*(J1-1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1) - &
246:   &               PARAM1*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1))/PARAM1))**2+ &246:   &               PARAM1*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+1)-XYZ(NOPT*J1+3*(J2-1)+1))/PARAM1))**2+ &
247:   &                           (XYZ(NOPT*(J1-1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2) - &247:   &                           (XYZ(NOPT*(J1-1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2) - &
248:   &               PARAM2*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2))/PARAM2))**2248:   &               PARAM2*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+2)-XYZ(NOPT*J1+3*(J2-1)+2))/PARAM2))**2
249:    IF (.NOT.TWOD) DISTB=DISTB+(XYZ(NOPT*(J1-1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3) - &249:    IF (.NOT.TWOD) DISTA=DISTA+(XYZ(NOPT*(J1-1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3) - &
250:   &               PARAM3*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3))/PARAM3))**2250:   &               PARAM3*NINT((XYZ(NOPT*(J1-1)+3*(J2-1)+3)-XYZ(NOPT*J1+3*(J2-1)+3))/PARAM3))**2
251:                      ENDDO251:                      ENDDO
252:                      SSS(NOPT*J1+1:NOPT*(J1+1))= - ( NEWNEBK(J1+1)*SQRT(DISTA) - NEWNEBK(J1)*SQRT(DISTB) )*TANVEC(1:NOPT,J1)252:                      SSS(NOPT*J1+1:NOPT*(J1+1))= - ( NEWNEBK(J1+1)*SQRT(DISTA) - NEWNEBK(J1)*SQRT(DISTB) )*TANVEC(1:NOPT,J1)
253:                   ELSE253:                   ELSE
254:                      SSS(NOPT*J1+1:NOPT*(J1+1))= &254:                      SSS(NOPT*J1+1:NOPT*(J1+1))= &
255:   &               - ( NEWNEBK(J1+1)*SQRT(SUM( ( XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) &255:   &               - ( NEWNEBK(J1+1)*SQRT(SUM( ( XYZ(NOPT*(J1+1)+1:NOPT*(J1+2)) - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) &
256:   &                 - NEWNEBK(J1)  *SQRT(SUM( ( XYZ(NOPT*(J1-1)+1:NOPT*J1)     - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) )*TANVEC(:,J1)256:   &                 - NEWNEBK(J1)  *SQRT(SUM( ( XYZ(NOPT*(J1-1)+1:NOPT*J1)     - XYZ(NOPT*J1+1:NOPT*(J1+1)) )**2 )) )*TANVEC(:,J1)
257:                   ENDIF257:                   ENDIF
258:                ENDDO258:                ENDDO
259:               CASE("dneb2") ! same as "dneb", except spring gradient usage is more consistent;259:               CASE("dneb2") ! same as "dneb", except spring gradient usage is more consistent;


r30520/key.f90 2016-05-26 15:30:14.304126387 +0100 r30519/key.f90 2016-05-26 15:30:16.528157411 +0100
 15:      &        NRBTRIES, REDOTSIM, REDOBFGSSTEPS, RPIMAGES, RPDOF, SDOXYGEN, SDHYDROGEN, SDCHARGE, BOWMANPES, & 15:      &        NRBTRIES, REDOTSIM, REDOBFGSSTEPS, RPIMAGES, RPDOF, SDOXYGEN, SDHYDROGEN, SDCHARGE, BOWMANPES, &
 16:      &        INTCONSEP, PATOM1, PATOM2, INTREPSEP, NCONSTRAINTON, CPREPSEP, CPCONSEP, NCONGEOM, & 16:      &        INTCONSEP, PATOM1, PATOM2, INTREPSEP, NCONSTRAINTON, CPREPSEP, CPCONSEP, NCONGEOM, &
 17:      &        NCPREPULSIVE, NCPCONSTRAINT, MAXCONUSE, INTCONSTEPS, INTRELSTEPS, INTSTEPS1, INTLJSTEPS, & 17:      &        NCPREPULSIVE, NCPCONSTRAINT, MAXCONUSE, INTCONSTEPS, INTRELSTEPS, INTSTEPS1, INTLJSTEPS, &
 18:      &        NTRAPPOW, MAXINTIMAGE, CHECKDID, CHECKREPINTERVAL, INTFREEZEMIN, INTNTRIESMAX, INTIMAGEINCR, & 18:      &        NTRAPPOW, MAXINTIMAGE, CHECKDID, CHECKREPINTERVAL, INTFREEZEMIN, INTNTRIESMAX, INTIMAGEINCR, &
 19:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, NRANROT, NENDDUP, LOCALPERMNEIGH, & 19:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, NRANROT, NENDDUP, LOCALPERMNEIGH, &
 20:      &        LOCALPERMMAXSEP, NONEDAPBC, STRUC, QCHEMESNAO, QCHEMESNMO, QCHEMESNZERO, QCHEMESNELEC, PMPATHINR, & 20:      &        LOCALPERMMAXSEP, NONEDAPBC, STRUC, QCHEMESNAO, QCHEMESNMO, QCHEMESNZERO, QCHEMESNELEC, PMPATHINR, &
 21:      &        MULTISUNIT, MULTIFUNIT,NIMAGEINST,NGLJ,ST_TSSTEP,LANSTEP,NONFREEZE, & 21:      &        MULTISUNIT, MULTIFUNIT,NIMAGEINST,NGLJ,ST_TSSTEP,LANSTEP,NONFREEZE, &
 22:      &        MCPATHBINS,MCPATHEQUIL,MCPATHSTEPS,MCPATHPRTFRQ,MCPATHTS,MCPATHSCHECK,RPHSLICES,RPHQBINS, & 22:      &        MCPATHBINS,MCPATHEQUIL,MCPATHSTEPS,MCPATHPRTFRQ,MCPATHTS,MCPATHSCHECK,RPHSLICES,RPHQBINS, &
 23:      &        ITWIST, JTWIST, KTWIST, LTWIST, MCPATHSTART, MCPATHBLOCK, MCPATHOVER, NCPU, MCPATHDOBLOCK, MCMERGES, MCMERGEQ, & 23:      &        ITWIST, JTWIST, KTWIST, LTWIST, MCPATHSTART, MCPATHBLOCK, MCPATHOVER, NCPU, MCPATHDOBLOCK, MCMERGES, MCMERGEQ, &
 24:      &        MCMERGEI,GAUSSIANCHARGE,GAUSSIANMULTI,ITG03, REDOTS, QCIPERMCHECKINT, & 24:      &        MCMERGEI,GAUSSIANCHARGE,GAUSSIANMULTI,ITG03, REDOTS, QCIPERMCHECKINT, &
 25:      &        MLPIN, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, N_TO_ALIGN, DJWRBID, STM, NHEXAMERS 25:      &        MLPIN, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, N_TO_ALIGN, DJWRBID, NHEXAMERS
 26:  26: 
 27:       LOGICAL :: DTEST, MASST, RTEST, EFSTEPST, VECTORST, SUMMARYT, DUMPV, DUMPMAG, FREEZE, FREEZERANGE, GRADSQ, & 27:       LOGICAL :: DTEST, MASST, RTEST, EFSTEPST, VECTORST, SUMMARYT, DUMPV, DUMPMAG, FREEZE, FREEZERANGE, GRADSQ, &
 28:      &        PGRAD, VALUEST, ADMT, BFGSMINT, BFGSTST, CHECKINDEX, TOSI, CONTAINER, & 28:      &        PGRAD, VALUEST, ADMT, BFGSMINT, BFGSTST, CHECKINDEX, TOSI, CONTAINER, &
 29:      &        GAUSSIAN, CADPAC, PRESSURE, FTEST, DCHECK, CP2K, DFTP, CPMD, CPMDC, FREEZERES, DF1T, & 29:      &        GAUSSIAN, CADPAC, PRESSURE, FTEST, DCHECK, CP2K, DFTP, CPMD, CPMDC, FREEZERES, DF1T, &
 30:      &        VARIABLES, FIELDT, OHT, IHT, TDT, D5HT, TWOENDS, PV, FRACTIONAL, BLNT, HYBRIDMINT, & 30:      &        VARIABLES, FIELDT, OHT, IHT, TDT, D5HT, TWOENDS, PV, FRACTIONAL, BLNT, HYBRIDMINT, &
 31:      &        INDEXT, LANCZOST, NOSHIFT, GAMESSUS, GAMESSUK, PVTS, RIGIDBODY, CASTEP, ONETEP, QCHEM, QCHEMES, VASP, & 31:      &        INDEXT, LANCZOST, NOSHIFT, GAMESSUS, GAMESSUK, PVTS, RIGIDBODY, CASTEP, ONETEP, QCHEM, QCHEMES, VASP, &
 32:      &        BFGSSTEP, BULKT, HUPDATE, NOHESS, READV, NOIT, THOMSONT, SIO2T, SIO2C6T, BISECTT, BISECTDEBUG, & 32:      &        BFGSSTEP, BULKT, HUPDATE, NOHESS, READV, NOIT, THOMSONT, SIO2T, SIO2C6T, BISECTT, BISECTDEBUG, &
 33:      &        TOSIC6, TOSIPOL, FIXIMAGE, DFTBT, CHECKCONT, CHECKDT, SHIFTED, READSP, DUMPSP, NOFRQS, & 33:      &        TOSIC6, TOSIPOL, FIXIMAGE, DFTBT, CHECKCONT, CHECKDT, SHIFTED, READSP, DUMPSP, NOFRQS, &
 34:      &        ALLSTEPS, ALLVECTORS, MWVECTORS, WELCH, BINARY, READHESS, MOVIE, NORESET, TWOD, & 34:      &        ALLSTEPS, ALLVECTORS, MWVECTORS, WELCH, BINARY, READHESS, MOVIE, NORESET, TWOD, &
 35:      &        DOUBLET, REOPT, PARALLEL, LINEMIN, FIXD, KEEPINDEX, BSMIN, PRINTPTS, RKMIN, REPELTST,& 35:      &        DOUBLET, REOPT, PARALLEL, LINEMIN, FIXD, KEEPINDEX, BSMIN, PRINTPTS, RKMIN, REPELTST,&
 45:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, & 45:      &        QSPCFWT, QTIP4PFT, CFUSIONT, DUMPINTXYZ, DUMPINTEOS, INTLJT, INTTST, EYTRAPT, OHCELLT, MKTRAPT, &
 46:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, & 46:      &        INTFREEZET, LPERMDIST, CHECKNEGATIVET, CHECKOVERLAPT, ACK1, ACK2, CONDATT, USERPOTT, &
 47:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, & 47:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, &
 48:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, & 48:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, &
 49:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 49:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 50:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 50:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 51:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 51:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 52:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 52:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 53:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 53:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 54:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 54:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 55:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, REBOXT 55:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT
 56:  56: 
 57: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 57: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 58:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 58:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 59:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 59:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 60:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 60:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
 61:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 61:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
 62:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads 62:       DOUBLE PRECISION :: HARMPOLY_BONLEN = 0.0d0 ! equilibrium length of springs between beads
 63:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs 63:       DOUBLE PRECISION :: HARMPOLY_K = 1.0d2 ! force constant of the springs
 64:  64: 
 65: ! hk286 > generalised THOMSON problem 65: ! hk286 > generalised THOMSON problem


r30520/keywords.f 2016-05-26 15:30:14.540129679 +0100 r30519/keywords.f 2016-05-26 15:30:16.760160649 +0100
898: 898: 
899:          CUDAT=.FALSE.899:          CUDAT=.FALSE.
900:          CUDAPOT=' '900:          CUDAPOT=' '
901:          CUDATIMET=.FALSE.901:          CUDATIMET=.FALSE.
902: 902: 
903:          CLIMBERT=.FALSE.903:          CLIMBERT=.FALSE.
904:          CLIMBERINIT=.FALSE.904:          CLIMBERINIT=.FALSE.
905:          CLIMBERSTEPS=20905:          CLIMBERSTEPS=20
906:          CLIMBERCONV=0.2906:          CLIMBERCONV=0.2
907:          CLIMBERSPRING=5.0907:          CLIMBERSPRING=5.0
908: ! 
909: ! Stealthy potential 
910: ! 
911:          STEALTHYT=.FALSE. 
912:          STEALTV=.FALSE. 
913: 908: 
914: !909: !
915: ! Neural network potential910: ! Neural network potential
916: !911: !
917:          MLP3T=.FALSE.912:          MLP3T=.FALSE.
918:          MLPB3T=.FALSE.913:          MLPB3T=.FALSE.
919:          MLPNEWREG=.FALSE.914:          MLPNEWREG=.FALSE.
920:          MLPPROB=.FALSE.915:          MLPPROB=.FALSE.
921:          MLPDONE=.FALSE.916:          MLPDONE=.FALSE.
922:          MLPNORM=.FALSE.917:          MLPNORM=.FALSE.
1613:                CALL READF(bulk_boxvec(2))1608:                CALL READF(bulk_boxvec(2))
1614:             ENDIF1609:             ENDIF
1615:             IF (NITEMS.GT.3) THEN1610:             IF (NITEMS.GT.3) THEN
1616:                CALL READF(bulk_boxvec(3))1611:                CALL READF(bulk_boxvec(3))
1617:             ENDIF1612:             ENDIF
1618: ! 1613: ! 
1619: ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC1614: ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1620: ! CADPAC tells the program to read derivative information in1615: ! CADPAC tells the program to read derivative information in
1621: ! CADPAC format.                                        - default FALSE1616: ! CADPAC format.                                        - default FALSE
1622: ! 1617: ! 
1623:          ELSE IF (WORD.EQ.'REBOX') THEN 
1624:             REBOXT=.TRUE. 
1625:  
1626:          ELSE IF (WORD.EQ.'CADPAC') THEN1618:          ELSE IF (WORD.EQ.'CADPAC') THEN
1627:             CADPAC=.TRUE.1619:             CADPAC=.TRUE.
1628:             CALL READA(SYS)1620:             CALL READA(SYS)
1629:             DO J1=1,801621:             DO J1=1,80
1630:                IF (SYS(J1:J1).EQ.' ') THEN1622:                IF (SYS(J1:J1).EQ.' ') THEN
1631:                   LSYS=J1-11623:                   LSYS=J1-1
1632:                   GOTO 101624:                   GOTO 10
1633:                ENDIF1625:                ENDIF
1634:             ENDDO1626:             ENDDO
1635: 10          IF (NITEMS.GT.2) THEN1627: 10          IF (NITEMS.GT.2) THEN
5653:             SSHT=.TRUE.5645:             SSHT=.TRUE.
5654: 5646: 
5655:          ELSE IF (WORD.EQ.'STEALTHY') THEN5647:          ELSE IF (WORD.EQ.'STEALTHY') THEN
5656:             STEALTHYT=.TRUE.5648:             STEALTHYT=.TRUE.
5657:             CALL READF(KLIM)5649:             CALL READF(KLIM)
5658:             IF (NITEMS.GT.2) THEN5650:             IF (NITEMS.GT.2) THEN
5659:                CALL READF(SCA)5651:                CALL READF(SCA)
5660:             ELSE5652:             ELSE
5661:                SCA=15653:                SCA=1
5662:             END IF5654:             END IF
5663:  
5664:          ELSE IF (WORD.EQ.'STTEST') THEN 
5665:             STEALTV=.TRUE. 
5666: ! 5655: ! 
5667: ! NSTEPMIN sets the minimum number of steps allowed before convergence.5656: ! NSTEPMIN sets the minimum number of steps allowed before convergence.
5668: ! 5657: ! 
5669:          ELSE IF (WORD .EQ. 'STEPMIN') THEN5658:          ELSE IF (WORD .EQ. 'STEPMIN') THEN
5670:             CALL READI(NSTEPMIN)5659:             CALL READI(NSTEPMIN)
5671: ! 5660: ! 
5672: ! STEPS n sets the number of optimisation steps to perform5661: ! STEPS n sets the number of optimisation steps to perform
5673: ! per call to OPTIM                                    - default n=15662: ! per call to OPTIM                                    - default n=1
5674: ! If BFGSSTEPS is not specified then it is set to the same value as NSTEPS5663: ! If BFGSSTEPS is not specified then it is set to the same value as NSTEPS
5675: ! 5664: ! 


r30520/lbfgs.f90 2016-05-26 15:30:13.644117180 +0100 r30519/lbfgs.f90 2016-05-26 15:30:15.872148263 +0100
 28:      USE NEBUTILS 28:      USE NEBUTILS
 29:      USE GRADIENTS 29:      USE GRADIENTS
 30:      USE NEBOUTPUT 30:      USE NEBOUTPUT
 31:      USE PORFUNCS 31:      USE PORFUNCS
 32:      USE MODCHARMM, ONLY : CHRMMT,CHECKOMEGAT 32:      USE MODCHARMM, ONLY : CHRMMT,CHECKOMEGAT
 33:      USE KEY, ONLY : FREEZENODEST, FREEZETOL, DESMDEBUG, MAXNEBBFGS, NEBMUPDATE, NEBDGUESS, & 33:      USE KEY, ONLY : FREEZENODEST, FREEZETOL, DESMDEBUG, MAXNEBBFGS, NEBMUPDATE, NEBDGUESS, &
 34:           & DESMAXEJUMP,INTEPSILON, DESMAXAVGE, KADJUSTFRAC, KADJUSTFRQ, DNEBSWITCH, KADJUSTTOL, NEBRESEEDT, & 34:           & DESMAXEJUMP,INTEPSILON, DESMAXAVGE, KADJUSTFRAC, KADJUSTFRQ, DNEBSWITCH, KADJUSTTOL, NEBRESEEDT, &
 35:           & NEBRESEEDINT, NEBRESEEDEMAX, NEBRESEEDBMAX, NEBKFINAL, NEBFACTOR, & 35:           & NEBRESEEDINT, NEBRESEEDEMAX, NEBRESEEDBMAX, NEBKFINAL, NEBFACTOR, &
 36:           & NREPMAX, ORDERI, ORDERJ, EPSALPHA, NREPULSIVE, DISTREF, ADDREPT, NEBKINITIAL, & 36:           & NREPMAX, ORDERI, ORDERJ, EPSALPHA, NREPULSIVE, DISTREF, ADDREPT, NEBKINITIAL, &
 37:           & NEBRESEEDDEL1, NEBRESEEDDEL2, NEBRESEEDPOW1, NEBRESEEDPOW2, REPPOW, & 37:           & NEBRESEEDDEL1, NEBRESEEDDEL2, NEBRESEEDPOW1, NEBRESEEDPOW2, REPPOW, &
 38:           & BULKT, REDOTSIM, NREPMAX, COLDFUSIONLIMIT, DNEBEFRAC,VARIABLES, MAXERISE 38:           & BULKT, REDOTSIM, NREPMAX, COLDFUSIONLIMIT, DNEBEFRAC, VARIABLES
 39:      USE INTCOMMONS, ONLY : DESMINT, NNZ, KD, NINTC, PREVDIH, DIHINFO 39:      USE INTCOMMONS, ONLY : DESMINT, NNZ, KD, NINTC, PREVDIH, DIHINFO
 40:      USE COMMONS, ONLY: REDOPATHNEB 40:      USE COMMONS, ONLY: REDOPATHNEB
 41: ! hk286 41: ! hk286
 42:      USE GENRIGID 42:      USE GENRIGID
 43:  43: 
 44:      IMPLICIT NONE  44:      IMPLICIT NONE 
 45:  45: 
 46:      INTEGER,INTENT(IN) :: D,U  ! DIMENSIONALITY OF THE PROBLEM AND NUMBER OF UPDATES 46:      INTEGER,INTENT(IN) :: D,U  ! DIMENSIONALITY OF THE PROBLEM AND NUMBER OF UPDATES
 47:      INTEGER NPERSIST           ! number of persistent minima 47:      INTEGER NPERSIST           ! number of persistent minima
 48:      INTEGER PERSISTTHRESH      ! persistence threshold 48:      INTEGER PERSISTTHRESH      ! persistence threshold
 49:      LOGICAL PERSISTENT(NIMAGE+2) ! logical array to identify persistent minima 49:      LOGICAL PERSISTENT(NIMAGE+2) ! logical array to identify persistent minima
 50:      INTEGER NTIMESMIN(NIMAGE+2)  ! number of consecutive steps this image is a local minimum 50:      INTEGER NTIMESMIN(NIMAGE+2)  ! number of consecutive steps this image is a local minimum
 51:      LOGICAL NOTNEW  51:      LOGICAL NOTNEW 
 52:      DOUBLE PRECISION DIJ, DMAX, DS, DF, EWORST, EMAX, BESTE, ENOLD 52:      DOUBLE PRECISION DIJ, DMAX, DS, DF, EWORST, EMAX, BESTE
 53:      DOUBLE PRECISION, ALLOCATABLE :: REPTEMP(:) 53:      DOUBLE PRECISION, ALLOCATABLE :: REPTEMP(:)
 54:      INTEGER, ALLOCATABLE :: IREPTEMP(:) 54:      INTEGER, ALLOCATABLE :: IREPTEMP(:)
 55:      INTEGER MAXIM, NDECREASE, NFAIL 55:      INTEGER MAXIM, NDECREASE, NFAIL
 56:      LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM 56:      LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM
 57:      COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 57:      COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 58:  58: 
 59:      INTEGER :: J2,POINT,BOUND,NPT,CP,I,ISTAT,J1,JMINUS,JPLUS,J3,JDO,J4,NPEPFAIL,NBADTOTAL 59:      INTEGER :: J2,POINT,BOUND,NPT,CP,I,ISTAT,J1,JMINUS,JPLUS,J3,JDO,J4,NPEPFAIL,NBADTOTAL
 60:      INTEGER AT1(NATOMS),AT2(NATOMS),AT3(NATOMS),AT4(NATOMS) 60:      INTEGER AT1(NATOMS),AT2(NATOMS),AT3(NATOMS),AT4(NATOMS)
 61:      DOUBLE PRECISION :: YS,YY,SQ,YR,BETA,GNORM,DDOT,DUMMY,STPMIN,PREVGRAD,MEANSEP,EMINUS,EPLUS, & 61:      DOUBLE PRECISION :: YS,YY,SQ,YR,BETA,GNORM,DDOT,DUMMY,STPMIN,PREVGRAD,MEANSEP,EMINUS,EPLUS, &
 62:   &                      LCOORDS(NOPT), ENERGY, INVDTOACTIVE(NATOMS), STIME 62:   &                      LCOORDS(NOPT), ENERGY, INVDTOACTIVE(NATOMS), STIME
 92:      IF (DUMPNEBXYZ) CALL RWG("w",.False.,0) 92:      IF (DUMPNEBXYZ) CALL RWG("w",.False.,0)
 93:      IF (DUMPNEBEOS) CALL WRITEPROFILE(0) 93:      IF (DUMPNEBEOS) CALL WRITEPROFILE(0)
 94:      IF (MOREPRINTING) THEN 94:      IF (MOREPRINTING) THEN
 95:           WRITE(*,'(A)') '                  Iter   Energy per image      RMS Force       Av.Dev    Path     Step' 95:           WRITE(*,'(A)') '                  Iter   Energy per image      RMS Force       Av.Dev    Path     Step'
 96:      ENDIF 96:      ENDIF
 97:      IF (PREVGRAD.LT.DNEBSWITCH) THEN 97:      IF (PREVGRAD.LT.DNEBSWITCH) THEN
 98:         CALL OLDNEBGRADIENT 98:         CALL OLDNEBGRADIENT
 99:      ELSE 99:      ELSE
100:         CALL NEBGRADIENT100:         CALL NEBGRADIENT
101:      ENDIF101:      ENDIF
102:  
103:      ENOLD=ETOTAL 
104:      IF (ETOTAL/NIMAGE.LT.COLDFUSIONLIMIT) THEN102:      IF (ETOTAL/NIMAGE.LT.COLDFUSIONLIMIT) THEN
105:         WRITE(*,'(A,2G20.10)') ' lbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/NIMAGE,COLDFUSIONLIMIT103:         WRITE(*,'(A,2G20.10)') ' lbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/NIMAGE,COLDFUSIONLIMIT
106:         IF (DEBUG) CALL DUMPFILES("e")104:         IF (DEBUG) CALL DUMPFILES("e")
107:         IF (KADJUSTFRQ.GT.0) PRINT '(A,G20.10)',' lbfgs> Final DNEB force constant ',NEWNEBK(1)105:         IF (KADJUSTFRQ.GT.0) PRINT '(A,G20.10)',' lbfgs> Final DNEB force constant ',NEWNEBK(1)
108:         RETURN106:         RETURN
109:      ENDIF107:      ENDIF
110: 108: 
111:      NITERDONE=1109:      NITERDONE=1
112: !    DO NITERDONE=1,MAX(NITERMAX,NITERMIN)110: !    DO NITERDONE=1,MAX(NITERMAX,NITERMIN)
113:      DO ! main do loop with counter NITERDONE111:      DO ! main do loop with counter NITERDONE
118:           IF (.NOT.DIAGCO) DIAG(1:D)=NEBDGUESS116:           IF (.NOT.DIAGCO) DIAG(1:D)=NEBDGUESS
119:           SEARCHSTEP(0,:)= -G*DIAG                           ! NR STEP FOR DIAGONAL INVERSE HESSIAN117:           SEARCHSTEP(0,:)= -G*DIAG                           ! NR STEP FOR DIAGONAL INVERSE HESSIAN
120:           GTMP           = SEARCHSTEP(0,:)118:           GTMP           = SEARCHSTEP(0,:)
121:           GNORM          = MAX(SQRT(DOT_PRODUCT(G,G)),1.0D-100)119:           GNORM          = MAX(SQRT(DOT_PRODUCT(G,G)),1.0D-100)
122:           STP            = MIN(1.0D0/GNORM, GNORM)           ! MAKE THE FIRST GUESS FOR THE STEP LENGTH CAUTIOUS120:           STP            = MIN(1.0D0/GNORM, GNORM)           ! MAKE THE FIRST GUESS FOR THE STEP LENGTH CAUTIOUS
123:      ELSE MAIN121:      ELSE MAIN
124:           BOUND=NITERDONE-1122:           BOUND=NITERDONE-1
125:           IF (NITERDONE.GT.NEBMUPDATE) BOUND=NEBMUPDATE123:           IF (NITERDONE.GT.NEBMUPDATE) BOUND=NEBMUPDATE
126:           YS=DOT_PRODUCT( GDIF(NPT/D,:), SEARCHSTEP(NPT/D,:)  )124:           YS=DOT_PRODUCT( GDIF(NPT/D,:), SEARCHSTEP(NPT/D,:)  )
127:           IF (YS==0.0D0) YS=1.0D0125:           IF (YS==0.0D0) YS=1.0D0
128:  126:       
129:           ! Update estimate of diagonal inverse Hessian elements.127:           ! Update estimate of diagonal inverse Hessian elements.
130:           ! We divide by both YS and YY at different points, so they had better not be zero!128:           ! We divide by both YS and YY at different points, so they had better not be zero!
131:           IF (.NOT.DIAGCO) THEN129:           IF (.NOT.DIAGCO) THEN
132:                YY=DOT_PRODUCT( GDIF(NPT/D,:) , GDIF(NPT/D,:) )130:                YY=DOT_PRODUCT( GDIF(NPT/D,:) , GDIF(NPT/D,:) )
133:                IF (YY==0.0D0) YY=1.0D0131:                IF (YY==0.0D0) YY=1.0D0
134: !              DIAG = ABS(YS/YY)132: !              DIAG = ABS(YS/YY)
135:                DIAG = YS/YY133:                DIAG = YS/YY
136:           ELSE134:           ELSE
137:                CALL CHECKINPUT135:                CALL CHECKINPUT
138:           ENDIF136:           ENDIF
139:  137:       
140:           ! COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, "Updating quasi-Newton matrices with limited storage",138:           ! COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, "Updating quasi-Newton matrices with limited storage",
141:           ! Mathematics of Computation, Vol.35, No.151, pp. 773-782139:           ! Mathematics of Computation, Vol.35, No.151, pp. 773-782
142:           CP= POINT; IF (POINT==0) CP = NEBMUPDATE140:           CP= POINT; IF (POINT==0) CP = NEBMUPDATE
143:           RHO1(CP)=1.0D0/YS141:           RHO1(CP)=1.0D0/YS
144:           GTMP = -G142:           GTMP = -G
145:           CP= POINT 143:           CP= POINT 
146:  144:                    
147:           DO I= 1,BOUND 145:           DO I= 1,BOUND 
148: !              CP = CP - 1; IF (CP == -1) CP = M - 1146: !              CP = CP - 1; IF (CP == -1) CP = M - 1
149:                CP = CP - 1; IF (CP == -1) CP = NEBMUPDATE - 1147:                CP = CP - 1; IF (CP == -1) CP = NEBMUPDATE - 1
150:                SQ= DOT_PRODUCT( SEARCHSTEP(CP,:),GTMP )148:                SQ= DOT_PRODUCT( SEARCHSTEP(CP,:),GTMP )
151:                ALPHA(CP+1) = RHO1(CP+1) * SQ149:                ALPHA(CP+1) = RHO1(CP+1) * SQ
152:                GTMP        = -ALPHA(CP+1)*GDIF(CP,:) + GTMP150:                GTMP        = -ALPHA(CP+1)*GDIF(CP,:) + GTMP
153:           ENDDO151:           ENDDO
154:               152:               
155:           GTMP=DIAG*GTMP153:           GTMP=DIAG*GTMP
156: 154: 
193:      DO J2=1,NIMAGE191:      DO J2=1,NIMAGE
194:           STEPIMAGE(J2) = SQRT(DOT_PRODUCT(SEARCHSTEP(POINT,NOPT*(J2-1)+1:NOPT*J2),SEARCHSTEP(POINT,NOPT*(J2-1)+1:NOPT*J2)))192:           STEPIMAGE(J2) = SQRT(DOT_PRODUCT(SEARCHSTEP(POINT,NOPT*(J2-1)+1:NOPT*J2),SEARCHSTEP(POINT,NOPT*(J2-1)+1:NOPT*J2)))
195:           DUMMY=STEPIMAGE(J2)193:           DUMMY=STEPIMAGE(J2)
196:           IF (STEPIMAGE(J2) > MAXNEBBFGS) THEN194:           IF (STEPIMAGE(J2) > MAXNEBBFGS) THEN
197:                STP(NOPT*(J2-1)+1:NOPT*J2) = MAXNEBBFGS/STEPIMAGE(J2)195:                STP(NOPT*(J2-1)+1:NOPT*J2) = MAXNEBBFGS/STEPIMAGE(J2)
198:                STPMIN=MIN(STPMIN,STP(NOPT*(J2-1)+1))196:                STPMIN=MIN(STPMIN,STP(NOPT*(J2-1)+1))
199:           ENDIF197:           ENDIF
200: !         PRINT '(A,I8,3G20.10)','image,initial step size,STP,prod=',J2,DUMMY,STP(NOPT*(J2-1)+1),STEPIMAGE(J2)*STP(NOPT*(J2-1)+1)198: !         PRINT '(A,I8,3G20.10)','image,initial step size,STP,prod=',J2,DUMMY,STP(NOPT*(J2-1)+1),STEPIMAGE(J2)*STP(NOPT*(J2-1)+1)
201:      ENDDO199:      ENDDO
202:      STP(1:D)=STPMIN200:      STP(1:D)=STPMIN
 201: 
203: ! EFK: decide whether to freeze some nodes202: ! EFK: decide whether to freeze some nodes
204:      IF (FREEZENODEST.AND.(.NOT.NEBRESEEDT)) THEN203:      IF (FREEZENODEST.AND.(.NOT.NEBRESEEDT)) THEN
205:         TOTGNORM = SQRT(DOT_PRODUCT(G(1:NOPT*NIMAGE),G(1:NOPT*NIMAGE))/NIMAGE)204:         TOTGNORM = SQRT(DOT_PRODUCT(G(1:NOPT*NIMAGE),G(1:NOPT*NIMAGE))/NIMAGE)
206:         DO IM = 1,NIMAGE205:         DO IM = 1,NIMAGE
207:            TESTG = SQRT(DOT_PRODUCT(G(NOPT*(IM-1)+1:NOPT*IM),G(NOPT*(IM-1)+1:NOPT*IM)))206:            TESTG = SQRT(DOT_PRODUCT(G(NOPT*(IM-1)+1:NOPT*IM),G(NOPT*(IM-1)+1:NOPT*IM)))
208:            IF(TESTG/TOTGNORM.LT.FREEZETOL) THEN207:            IF(TESTG/TOTGNORM.LT.FREEZETOL) THEN
209:               IF (MOREPRINTING) PRINT '(A,I6,2G20.10)', ' lbfgs> Freezing image: ', IM, TESTG, TOTGNORM208:               IF (MOREPRINTING) PRINT '(A,I6,2G20.10)', ' lbfgs> Freezing image: ', IM, TESTG, TOTGNORM
210:               IMGFREEZE(IM)=.TRUE.209:               IMGFREEZE(IM)=.TRUE.
211:               STEPIMAGE(IM)=0.0D0210:               STEPIMAGE(IM)=0.0D0
212:               STP(NOPT*(IM-1)+1:NOPT*IM)=0.0D0211:               STP(NOPT*(IM-1)+1:NOPT*IM)=0.0D0
256:         ENDDO255:         ENDDO
257:      ENDIF256:      ENDIF
258:      NDECREASE=0257:      NDECREASE=0
259: 20   X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)258: 20   X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)
260: 259: 
261:      IF (PREVGRAD.LT.DNEBSWITCH) THEN260:      IF (PREVGRAD.LT.DNEBSWITCH) THEN
262:         CALL OLDNEBGRADIENT261:         CALL OLDNEBGRADIENT
263:      ELSE262:      ELSE
264:         CALL NEBGRADIENT263:         CALL NEBGRADIENT
265:      ENDIF264:      ENDIF
266:  
267:      DO WHILE ((ETOTAL-ENOLD).GT.MAXERISE) 
268:         X(1:D) = X(1:D) - STP(1:D)*SEARCHSTEP(POINT,1:D) 
269:         STP(1:D) = STP(1:D)/10 
270:         X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D) 
271:         IF (PREVGRAD.LT.DNEBSWITCH) THEN 
272:            CALL OLDNEBGRADIENT 
273:         ELSE 
274:            CALL NEBGRADIENT 
275:         END IF 
276:      END DO 
277:  
278:      ENOLD=ETOTAL 
279:  
280:      IF (ETOTAL/NIMAGE.LT.COLDFUSIONLIMIT) THEN265:      IF (ETOTAL/NIMAGE.LT.COLDFUSIONLIMIT) THEN
281:         WRITE(*,'(A,2G20.10)') ' lbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/NIMAGE,COLDFUSIONLIMIT266:         WRITE(*,'(A,2G20.10)') ' lbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/NIMAGE,COLDFUSIONLIMIT
282:         IF (DEBUG) CALL DUMPFILES("e")267:         IF (DEBUG) CALL DUMPFILES("e")
283:         IF (KADJUSTFRQ.GT.0) PRINT '(A,G20.10)',' lbfgs> Final DNEB force constant ',NEWNEBK(1)268:         IF (KADJUSTFRQ.GT.0) PRINT '(A,G20.10)',' lbfgs> Final DNEB force constant ',NEWNEBK(1)
284:         RETURN269:         RETURN
285:      ENDIF270:      ENDIF
286: 271: 
287:      CALL IMAGEDISTRIBUTION272:      CALL IMAGEDISTRIBUTION
288: !273: !
289: !  Dynamic adjustment of NEWNEBK values. Local values are changed by KADJUSTFRAC if the274: !  Dynamic adjustment of NEWNEBK values. Local values are changed by KADJUSTFRAC if the
315:      ELSE300:      ELSE
316: !301: !
317: !  Alternative turning on of NEWNEBK.302: !  Alternative turning on of NEWNEBK.
318: !303: !
319:         DO J1=1,NIMAGE+1304:         DO J1=1,NIMAGE+1
320:            NEWNEBK(J1)=MIN(NEBKFINAL,NEWNEBK(J1)*NEBFACTOR)305:            NEWNEBK(J1)=MIN(NEBKFINAL,NEWNEBK(J1)*NEBFACTOR)
321:         ENDDO306:         ENDDO
322: !       IF (DEBUG) PRINT '(A,G20.10)','lbfgs> DNEB force constant for image 1 is ',NEWNEBK(1)307: !       IF (DEBUG) PRINT '(A,G20.10)','lbfgs> DNEB force constant for image 1 is ',NEWNEBK(1)
323:      ENDIF308:      ENDIF
324: 309: 
325:      STEPTOT = SUM(STEPIMAGE*STP)/NIMAGE310:      STEPTOT = SUM(STEPIMAGE)/NIMAGE
326: 311: 
327:      IF (MOREPRINTING) THEN312:      IF (MOREPRINTING) THEN
328:         IF (PREVGRAD.LT.DNEBSWITCH) THEN313:         IF (PREVGRAD.LT.DNEBSWITCH) THEN
329:            WRITE(*,'(A,I6,2G20.10,F8.2,A,F8.3,F9.3,A)') ' lbfgs> steps: ',NITERDONE,ETOTAL/NIMAGE,RMS,AVDEV,'%',SEPARATION, &314:            WRITE(*,'(A,I6,2G20.10,F8.2,A,F8.3,F9.3,A)') ' lbfgs> steps: ',NITERDONE,ETOTAL/NIMAGE,RMS,AVDEV,'%',SEPARATION, &
330:   &                                                     STEPTOT, ' singly-nudged'315:   &                                                     STEPTOT, ' singly-nudged'
331:         ELSE316:         ELSE
332:            WRITE(*,'(A,I6,2G20.10,F8.2,A,F8.3,F18.15)') ' lbfgs> steps: ',NITERDONE,ETOTAL/NIMAGE,RMS,AVDEV,'%',SEPARATION, STEPTOT317:            WRITE(*,'(A,I6,2G20.10,F8.2,A,F8.3,F9.3)') ' lbfgs> steps: ',NITERDONE,ETOTAL/NIMAGE,RMS,AVDEV,'%',SEPARATION, STEPTOT
333:         ENDIF318:         ENDIF
334:         CALL FLUSH(6)319:         CALL FLUSH(6)
335:      ENDIF320:      ENDIF
336:      IF (ETOTAL/NIMAGE.LT.BESTE) BESTE=ETOTAL/NIMAGE321:      IF (ETOTAL/NIMAGE.LT.BESTE) BESTE=ETOTAL/NIMAGE
337:      IF (DEBUG) CALL DUMPFILES("m")322:      IF (DEBUG) CALL DUMPFILES("m")
338: 323: 
339: 324: 
340:      CALL TERMINATERUN325:      CALL TERMINATERUN
341: 326: 
342: !327: !


r30520/minpermdist.f90 2016-05-26 15:30:14.760132749 +0100 r30519/minpermdist.f90 2016-05-26 15:30:16.972163606 +0100
 50: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the 50: !  The centres of coordinates for COORDSA and COORDSB can be anywhere. On return, the
 51: !  centre of coordinates of COORDSA will be the same as for COORDSB, unless we 51: !  centre of coordinates of COORDSA will be the same as for COORDSB, unless we
 52: !  are doing an ion trap potential. 52: !  are doing an ion trap potential.
 53: ! 53: !
 54: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST) 54: SUBROUTINE MINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,RIGID,RMATBEST)
 55: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, STOCKT, GEOMDIFFTOL, AMBERT, & 55: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, STOCKT, GEOMDIFFTOL, AMBERT, &
 56:   &            NFREEZE, NABT, RBAAT, ANGLEAXIS2, BESTPERM, LOCALPERMDIST, PULLT, EFIELDT, NTSITES, & 56:   &            NFREEZE, NABT, RBAAT, ANGLEAXIS2, BESTPERM, LOCALPERMDIST, PULLT, EFIELDT, NTSITES, &
 57:   &            RIGIDBODY, PERMDIST, OHCELLT, LPERMDIST, EYTRAPT, MKTRAPT, LOCALPERMCUT, LOCALPERMCUT2, & 57:   &            RIGIDBODY, PERMDIST, OHCELLT, LPERMDIST, EYTRAPT, MKTRAPT, LOCALPERMCUT, LOCALPERMCUT2, &
 58:   &            LOCALPERMCUTINC, NOINVERSION, MIEFT, & 58:   &            LOCALPERMCUTINC, NOINVERSION, MIEFT, &
 59:   &            EDIFFTOL, GMAX, CONVR, ATOMMATCHDIST, NRANROT, GTHOMSONT, GTHOMMET, & ! hk286 59:   &            EDIFFTOL, GMAX, CONVR, ATOMMATCHDIST, NRANROT, GTHOMSONT, GTHOMMET, & ! hk286
 60:   &            PHI4MODT, MCPATHT, AMBER12T, VARIABLES, MKTRAPT, ALIGNRBST, REBOXT 60:   &            PHI4MODT, MCPATHT, AMBER12T, VARIABLES, MKTRAPT, ALIGNRBST
 61: USE COMMONS,ONLY : NOPT 61: USE COMMONS,ONLY : NOPT
 62: USE MODCHARMM,ONLY : CHRMMT 62: USE MODCHARMM,ONLY : CHRMMT
 63: USE MODAMBER9, ONLY: NOPERMPROCHIRAL, PROCHIRALH 63: USE MODAMBER9, ONLY: NOPERMPROCHIRAL, PROCHIRALH
 64: USE INTCOMMONS, ONLY : INTMINPERMT, INTINTERPT, DESMINT, OLDINTMINPERMT, INTDISTANCET 64: USE INTCOMMONS, ONLY : INTMINPERMT, INTINTERPT, DESMINT, OLDINTMINPERMT, INTDISTANCET
 65: USE INTCUTILS, ONLY : INTMINPERM, OLD_INTMINPERM, INTMINPERM_CHIRAL, INTDISTANCE 65: USE INTCUTILS, ONLY : INTMINPERM, OLD_INTMINPERM, INTMINPERM_CHIRAL, INTDISTANCE
 66: USE GENRIGID 66: USE GENRIGID
 67: USE AMBER12_INTERFACE_MOD 67: USE AMBER12_INTERFACE_MOD
 68: USE CHIRALITY 68: USE CHIRALITY
 69: IMPLICIT NONE 69: IMPLICIT NONE
 70:  70: 
360:       CMAX=CMAX/NATOMS; CMAY=CMAY/NATOMS; CMAZ=CMAZ/NATOMS360:       CMAX=CMAX/NATOMS; CMAY=CMAY/NATOMS; CMAZ=CMAZ/NATOMS
361:       CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0361:       CMBX=0.0D0; CMBY=0.0D0; CMBZ=0.0D0
362:       DO J1=1,NATOMS362:       DO J1=1,NATOMS
363:          CMBX=CMBX+COORDSB(3*(J1-1)+1)363:          CMBX=CMBX+COORDSB(3*(J1-1)+1)
364:          CMBY=CMBY+COORDSB(3*(J1-1)+2)364:          CMBY=CMBY+COORDSB(3*(J1-1)+2)
365:          CMBZ=CMBZ+COORDSB(3*(J1-1)+3)365:          CMBZ=CMBZ+COORDSB(3*(J1-1)+3)
366:       ENDDO366:       ENDDO
367:       CMBX=CMBX/NATOMS; CMBY=CMBY/NATOMS; CMBZ=CMBZ/NATOMS367:       CMBX=CMBX/NATOMS; CMBY=CMBY/NATOMS; CMBZ=CMBZ/NATOMS
368:    ENDIF368:    ENDIF
369: ENDIF369: ENDIF
370:  
371: !370: !
372: ! It is possible for the standard orientation to result in a distance that is worse than371: ! It is possible for the standard orientation to result in a distance that is worse than
373: ! the starting distance. Hence we need to set XBEST here.372: ! the starting distance. Hence we need to set XBEST here.
374: !373: !
375: 374: 
376: DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS)375: DUMMYA(1:3*NATOMS)=COORDSA(1:3*NATOMS)
377: DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS)376: DUMMYB(1:3*NATOMS)=COORDSB(1:3*NATOMS)
378: IF (REBOXT) THEN 
379:    DO J1=1,NATOMS 
380:       DUMMYA(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-BOXLX*NINT(DUMMYA(3*(J1-1)+1)/BOXLX) 
381:       DUMMYA(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-BOXLY*NINT(DUMMYA(3*(J1-1)+2)/BOXLY) 
382:       IF (.NOT.TWOD) DUMMYA(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-BOXLZ*NINT(DUMMYA(3*(J1-1)+3)/BOXLZ) 
383:    ENDDO 
384:    DO J1=1,NATOMS 
385:       DUMMYB(3*(J1-1)+1)=DUMMYB(3*(J1-1)+1)-BOXLX*NINT(DUMMYB(3*(J1-1)+1)/BOXLX) 
386:       DUMMYB(3*(J1-1)+2)=DUMMYB(3*(J1-1)+2)-BOXLY*NINT(DUMMYB(3*(J1-1)+2)/BOXLY) 
387:       IF (.NOT.TWOD) DUMMYB(3*(J1-1)+3)=DUMMYB(3*(J1-1)+3)-BOXLZ*NINT(DUMMYB(3*(J1-1)+3)/BOXLZ) 
388:    ENDDO 
389: ENDIF 
390: DBEST=1.0D100377: DBEST=1.0D100
391: CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)378: CALL NEWMINDIST(DUMMYB,DUMMYA,NATOMS,DISTANCE,BULKT,TWOD,'AX    ',.FALSE.,RIGID,DEBUG,RMAT)
392: DBEST=DISTANCE**2379: DBEST=DISTANCE**2
393: 380: 
394: IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> Initial distance before standard orientation=',DISTANCE381: IF (DEBUG) PRINT '(A,G20.10)',' minpermdist> Initial distance before standard orientation=',DISTANCE
395: 382: 
396: ! If global translation is possible (i.e. no frozen atoms and not periodic boundary conditions) then383: ! If global translation is possible (i.e. no frozen atoms and not periodic boundary conditions) then
397: ! translate the origin of structure A to the CoM of structure B.384: ! translate the origin of structure A to the CoM of structure B.
398: IF ((NFREEZE.GT.0).OR.BULKT.OR.MIEFT.OR.MKTRAPT) THEN385: IF ((NFREEZE.GT.0).OR.BULKT.OR.MIEFT.OR.MKTRAPT) THEN
399:    XBEST(1:3*NATOMS)=DUMMYA(1:3*NATOMS)386:    XBEST(1:3*NATOMS)=DUMMYA(1:3*NATOMS)
403:       XBEST(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-CMBY390:       XBEST(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-CMBY
404:       XBEST(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-CMBZ391:       XBEST(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-CMBZ
405:    ENDDO392:    ENDDO
406: ELSE393: ELSE
407:    DO J1=1,NATOMS394:    DO J1=1,NATOMS
408:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX395:       XBEST(3*(J1-1)+1)=DUMMYA(3*(J1-1)+1)-CMBX
409:       XBEST(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-CMBY396:       XBEST(3*(J1-1)+2)=DUMMYA(3*(J1-1)+2)-CMBY
410:       XBEST(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-CMBZ397:       XBEST(3*(J1-1)+3)=DUMMYA(3*(J1-1)+3)-CMBZ
411:    ENDDO398:    ENDDO
412: ENDIF399: ENDIF
413:  
414: BESTINVERT=1400: BESTINVERT=1
415: DO J1=1,NATOMS401: DO J1=1,NATOMS
416:    BESTPERM(J1)=J1402:    BESTPERM(J1)=J1
417: ENDDO403: ENDDO
418: RMATBEST(1:3,1:3)=RMAT(1:3,1:3)404: RMATBEST(1:3,1:3)=RMAT(1:3,1:3)
419: ROTINVBBEST(1:3,1:3)=0.0D0405: ROTINVBBEST(1:3,1:3)=0.0D0
420: ROTINVBBEST(1,1)=1.0D0;ROTINVBBEST(2,2)=1.0D0;ROTINVBBEST(3,3)=1.0D0;406: ROTINVBBEST(1,1)=1.0D0;ROTINVBBEST(2,2)=1.0D0;ROTINVBBEST(3,3)=1.0D0;
421: ROTABEST(1:3,1:3)=0.0D0407: ROTABEST(1:3,1:3)=0.0D0
422: ROTABEST(1,1)=1.0D0;ROTABEST(2,2)=1.0D0;ROTABEST(3,3)=1.0D0;408: ROTABEST(1,1)=1.0D0;ROTABEST(2,2)=1.0D0;ROTABEST(3,3)=1.0D0;
423: !409: !
742:    ENDDO728:    ENDDO
743: ENDIF729: ENDIF
744: !730: !
745: ! Update coordinates in DUMMYA to overall permutation using NEWPERM.731: ! Update coordinates in DUMMYA to overall permutation using NEWPERM.
746: !732: !
747: DO J3=1,NATOMS733: DO J3=1,NATOMS
748:    DUMMYA(3*(J3-1)+1)=DUMMY(3*(NEWPERM(J3)-1)+1)734:    DUMMYA(3*(J3-1)+1)=DUMMY(3*(NEWPERM(J3)-1)+1)
749:    DUMMYA(3*(J3-1)+2)=DUMMY(3*(NEWPERM(J3)-1)+2)735:    DUMMYA(3*(J3-1)+2)=DUMMY(3*(NEWPERM(J3)-1)+2)
750:    DUMMYA(3*(J3-1)+3)=DUMMY(3*(NEWPERM(J3)-1)+3)736:    DUMMYA(3*(J3-1)+3)=DUMMY(3*(NEWPERM(J3)-1)+3)
751:    IF (J3.NE.NEWPERM(J3)) THEN737:    IF (J3.NE.NEWPERM(J3)) THEN
752:       !IF (DEBUG) WRITE(*,'(A,I5,A,I5)') ' minpermdist> move position ',NEWPERM(J3),' to ',J3738: !     IF (DEBUG) WRITE(*,'(A,I5,A,I5)') ' minpermdist> move position ',NEWPERM(J3),' to ',J3
753:       NPERM=NPERM+1739:       NPERM=NPERM+1
754:    ENDIF740:    ENDIF
755:    IF (STOCKT.OR.ANGLEAXIS2) THEN741:    IF (STOCKT.OR.ANGLEAXIS2) THEN
756:       IF (J3.LE.(NATOMS/2)) THEN742:       IF (J3.LE.(NATOMS/2)) THEN
757:          DISTANCE=DISTANCE+(DUMMYA(3*(J3-1)+1)-DUMMYB(3*(J3-1)+1))**2 &743:          DISTANCE=DISTANCE+(DUMMYA(3*(J3-1)+1)-DUMMYB(3*(J3-1)+1))**2 &
758:   &                       +(DUMMYA(3*(J3-1)+2)-DUMMYB(3*(J3-1)+2))**2 &744:   &                       +(DUMMYA(3*(J3-1)+2)-DUMMYB(3*(J3-1)+2))**2 &
759:   &                       +(DUMMYA(3*(J3-1)+3)-DUMMYB(3*(J3-1)+3))**2745:   &                       +(DUMMYA(3*(J3-1)+3)-DUMMYB(3*(J3-1)+3))**2
760:       ENDIF746:       ENDIF
761:    ELSEIF (.NOT.BULKT) THEN747:    ELSEIF (.NOT.BULKT) THEN
762:       DISTANCE=DISTANCE+(DUMMYA(3*(J3-1)+1)-DUMMYB(3*(J3-1)+1))**2 &748:       DISTANCE=DISTANCE+(DUMMYA(3*(J3-1)+1)-DUMMYB(3*(J3-1)+1))**2 &
872:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &858:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &
873:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &859:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
874:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2860:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
875:       ENDDO861:       ENDDO
876:    ELSE IF (STOCKT) THEN862:    ELSE IF (STOCKT) THEN
877:       CALL NEWROTGEOMSTOCK(NATOMS,XBEST,ROTINVBBEST,0.0D0,0.0D0,0.0D0)863:       CALL NEWROTGEOMSTOCK(NATOMS,XBEST,ROTINVBBEST,0.0D0,0.0D0,0.0D0)
878:       XDUMMY=0.0D0864:       XDUMMY=0.0D0
879:       DO J1=1,(NATOMS/2)865:       DO J1=1,(NATOMS/2)
880:          XBEST(3*(J1-1)+1)=XBEST(3*(J1-1)+1)+CMBX866:          XBEST(3*(J1-1)+1)=XBEST(3*(J1-1)+1)+CMBX
881:          XBEST(3*(J1-1)+2)=XBEST(3*(J1-1)+2)+CMBY867:          XBEST(3*(J1-1)+2)=XBEST(3*(J1-1)+2)+CMBY
882:          xbest(3*(J1-1)+3)=XBEST(3*(J1-1)+3)+CMBZ868:          XBEST(3*(J1-1)+3)=XBEST(3*(J1-1)+3)+CMBZ
883:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &869:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))**2+ &
884:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &870:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))**2+ &
885:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2871:   &                    (COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3))**2
886:       ENDDO872:       ENDDO
887:    ELSEIF (BULKT) THEN873:    ELSEIF (BULKT) THEN
888:       XDUMMY=0.0D0874:       XDUMMY=0.0D0
889:       DO J1=1,NATOMS875:       DO J1=1,NATOMS
890:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1) - BOXLX*NINT((COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))/BOXLX))**2+ &876:          XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1) - BOXLX*NINT((COORDSB(3*(J1-1)+1)-XBEST(3*(J1-1)+1))/BOXLX))**2+ &
891:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2) - BOXLY*NINT((COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))/BOXLY))**2877:   &                    (COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2) - BOXLY*NINT((COORDSB(3*(J1-1)+2)-XBEST(3*(J1-1)+2))/BOXLY))**2
892:          IF (.NOT.TWOD) XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3) - &878:          IF (.NOT.TWOD) XDUMMY=XDUMMY+(COORDSB(3*(J1-1)+3)-XBEST(3*(J1-1)+3) - &


r30520/newmindist.f90 2016-05-26 15:30:14.980135818 +0100 r30519/newmindist.f90 2016-05-26 15:30:17.208166899 +0100
 25: ! 25: !
 26: ! jmc As long as zsym isn't 'W' (in which case mind does something special) mind 26: ! jmc As long as zsym isn't 'W' (in which case mind does something special) mind
 27: ! doesn't care what atomic symbol we give it. 27: ! doesn't care what atomic symbol we give it.
 28: ! 28: !
 29: !  If PRESERVET is false we put RB into best correspondence with RA. This involves 29: !  If PRESERVET is false we put RB into best correspondence with RA. This involves
 30: !  a translation to the same centre of coordinates, followed by a rotation about that 30: !  a translation to the same centre of coordinates, followed by a rotation about that
 31: !  centre. 31: !  centre.
 32: ! 32: !
 33: SUBROUTINE NEWMINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET,RIGIDBODY,DEBUG,RMAT) 33: SUBROUTINE NEWMINDIST(RA,RB,NATOMS,DIST,BULKT,TWOD,ZUSE,PRESERVET,RIGIDBODY,DEBUG,RMAT)
 34: USE COMMONS,ONLY : PARAM1, PARAM2, PARAM3, NOPT 34: USE COMMONS,ONLY : PARAM1, PARAM2, PARAM3, NOPT
 35: USE KEY,ONLY : STOCKT, NFREEZE, MIEFT, RBAAT, PULLT, EFIELDT, EYTRAPT, MKTRAPT, USERPOTT, & 35: USE KEY,ONLY : STOCKT, NFREEZE, MIEFT, RBAAT, PULLT, EFIELDT, EYTRAPT, MKTRAPT, USERPOTT, GTHOMSONT, VARIABLES ! hk286
 36:   &            GTHOMSONT, VARIABLES, REBOXT ! hk286 
 37:  36: 
 38: IMPLICIT NONE 37: IMPLICIT NONE
 39:  38: 
 40: INTEGER J1, NATOMS, NSIZE, INFO, JINFO, JMIN 39: INTEGER J1, NATOMS, NSIZE, INFO, JINFO, JMIN
 41: DOUBLE PRECISION RA(NOPT), RB(NOPT), DIST, QMAT(4,4), XM, YM, ZM, XP, YP, ZP, OVEC(3), H1VEC(3), H2VEC(3), & 40: DOUBLE PRECISION RA(NOPT), RB(NOPT), DIST, QMAT(4,4), XM, YM, ZM, XP, YP, ZP, OVEC(3), H1VEC(3), H2VEC(3), &
 42:   &              DIAG(4), TEMPA(3*NOPT), RMAT(3,3), MINV, Q1, Q2, Q3, Q4, CMXA, CMYA, CMZA, CMXB, CMYB, CMZB, & 41:   &              DIAG(4), TEMPA(3*NOPT), RMAT(3,3), MINV, Q1, Q2, Q3, Q4, CMXA, CMYA, CMZA, CMXB, CMYB, CMZB, &
 43:   &              MYROTMAT(3,3), OMEGATOT(3,3) 42:   &              MYROTMAT(3,3), OMEGATOT(3,3)
 44: DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:) 43: DOUBLE PRECISION, ALLOCATABLE :: XA(:), XB(:)
 45: LOGICAL BULKT, TWOD, RIGIDBODY, PRESERVET, DEBUG, UPRETURN 44: LOGICAL BULKT, TWOD, RIGIDBODY, PRESERVET, DEBUG, UPRETURN
 46: CHARACTER(LEN=5) ZUSE 45: CHARACTER(LEN=5) ZUSE
162:    &               + (XA(3*(J1-1)+2)-XB(3*(J1-1)+2))**2 &161:    &               + (XA(3*(J1-1)+2)-XB(3*(J1-1)+2))**2 &
163:    &               + (XA(3*(J1-1)+3)-XB(3*(J1-1)+3))**2162:    &               + (XA(3*(J1-1)+3)-XB(3*(J1-1)+3))**2
164:       ENDDO163:       ENDDO
165:    ENDIF164:    ENDIF
166:    DIST=SQRT(DIST)165:    DIST=SQRT(DIST)
167:    166:    
168:    RMAT(1:3,1:3)=0.0D0 ! rotation matrix is the identity167:    RMAT(1:3,1:3)=0.0D0 ! rotation matrix is the identity
169:    RMAT(1,1)=1.0D0; RMAT(2,2)=1.0D0; RMAT(3,3)=1.0D0168:    RMAT(1,1)=1.0D0; RMAT(2,2)=1.0D0; RMAT(3,3)=1.0D0
170:    RETURN169:    RETURN
171: ENDIF170: ENDIF
172: IF (BULKT.AND. .NOT. REBOXT) THEN171: IF (BULKT) THEN
173:    BOXLX=PARAM1; BOXLY=PARAM2; BOXLZ=PARAM3172:    BOXLX=PARAM1; BOXLY=PARAM2; BOXLZ=PARAM3
174:    DO J1=1,NSIZE173:    DO J1=1,NSIZE
175:       XA(3*(J1-1)+1)=XA(3*(J1-1)+1)-BOXLX*NINT(XA(3*(J1-1)+1)/BOXLX)174:       XA(3*(J1-1)+1)=XA(3*(J1-1)+1)-BOXLX*NINT(XA(3*(J1-1)+1)/BOXLX)
176:       XA(3*(J1-1)+2)=XA(3*(J1-1)+2)-BOXLY*NINT(XA(3*(J1-1)+2)/BOXLY)175:       XA(3*(J1-1)+2)=XA(3*(J1-1)+2)-BOXLY*NINT(XA(3*(J1-1)+2)/BOXLY)
177:       IF (.NOT.TWOD) XA(3*(J1-1)+3)=XA(3*(J1-1)+3)-BOXLZ*NINT(XA(3*(J1-1)+3)/BOXLZ)176:       IF (.NOT.TWOD) XA(3*(J1-1)+3)=XA(3*(J1-1)+3)-BOXLZ*NINT(XA(3*(J1-1)+3)/BOXLZ)
178:    ENDDO177:    ENDDO
179:    DO J1=1,NSIZE178:    DO J1=1,NSIZE
180:       XB(3*(J1-1)+1)=XB(3*(J1-1)+1)-BOXLX*NINT(XB(3*(J1-1)+1)/BOXLX)179:       XB(3*(J1-1)+1)=XB(3*(J1-1)+1)-BOXLX*NINT(XB(3*(J1-1)+1)/BOXLX)
181:       XB(3*(J1-1)+2)=XB(3*(J1-1)+2)-BOXLY*NINT(XB(3*(J1-1)+2)/BOXLY)180:       XB(3*(J1-1)+2)=XB(3*(J1-1)+2)-BOXLY*NINT(XB(3*(J1-1)+2)/BOXLY)
182:       IF (.NOT.TWOD) XB(3*(J1-1)+3)=XB(3*(J1-1)+3)-BOXLZ*NINT(XB(3*(J1-1)+3)/BOXLZ)181:       IF (.NOT.TWOD) XB(3*(J1-1)+3)=XB(3*(J1-1)+3)-BOXLZ*NINT(XB(3*(J1-1)+3)/BOXLZ)
400:       CALL NEWROTGEOM(NSIZE,RB,RMAT,0.0D0,0.0D0,0.0D0)399:       CALL NEWROTGEOM(NSIZE,RB,RMAT,0.0D0,0.0D0,0.0D0)
401:    ELSE400:    ELSE
402: !401: !
403: !  Translate the RB coordinates to the centre of coordinates of RA.402: !  Translate the RB coordinates to the centre of coordinates of RA.
404: !403: !
405:       DO J1=1,NSIZE404:       DO J1=1,NSIZE
406:          RB(3*(J1-1)+1)=RB(3*(J1-1)+1)-CMXB+CMXA+XSHIFT405:          RB(3*(J1-1)+1)=RB(3*(J1-1)+1)-CMXB+CMXA+XSHIFT
407:          RB(3*(J1-1)+2)=RB(3*(J1-1)+2)-CMYB+CMYA+YSHIFT406:          RB(3*(J1-1)+2)=RB(3*(J1-1)+2)-CMYB+CMYA+YSHIFT
408:          RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-CMZB+CMZA+ZSHIFT407:          RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-CMZB+CMZA+ZSHIFT
409:       ENDDO408:       ENDDO
410:       !write(*,*) "movetest" 
411:       !write(*,*) -CMXB+CMXA+XSHIFT, -CMYB+CMYA+YSHIFT, -CMZB+CMZA+ZSHIFT 
412:  
413:       !DO J1=1,NSIZE 
414:       !   RB(3*(J1-1)+1)=RB(3*(J1-1)+1)-BOXLX*NINT(RB(3*(J1-1)+1)/BOXLX) 
415:       !   RB(3*(J1-1)+2)=RB(3*(J1-1)+2)-BOXLY*NINT(RB(3*(J1-1)+2)/BOXLY) 
416:       !   RB(3*(J1-1)+3)=RB(3*(J1-1)+3)-BOXLZ*NINT(RB(3*(J1-1)+3)/BOXLZ) 
417:       !ENDDO 
418: 409: 
419:       IF (.NOT.BULKT) THEN410:       IF (.NOT.BULKT) THEN
420:          IF (STOCKT) THEN411:          IF (STOCKT) THEN
421:             CALL NEWROTGEOMSTOCK(NATOMS,RB,RMAT,CMXA,CMYA,CMZA)412:             CALL NEWROTGEOMSTOCK(NATOMS,RB,RMAT,CMXA,CMYA,CMZA)
422:          ELSE413:          ELSE
423:             CALL NEWROTGEOM(NSIZE,RB,RMAT,CMXA,CMYA,CMZA)414:             CALL NEWROTGEOM(NSIZE,RB,RMAT,CMXA,CMYA,CMZA)
424:          ENDIF415:          ENDIF
425:       ENDIF416:       ENDIF
426:    ENDIF417:    ENDIF
427: ENDIF418: ENDIF


r30520/nnutils.f90 2016-05-26 15:30:13.860120193 +0100 r30519/nnutils.f90 2016-05-26 15:30:16.084151218 +0100
 26:           IMPLICIT NONE 26:           IMPLICIT NONE
 27:       27:      
 28:           AVDEVOLD=AVDEV 28:           AVDEVOLD=AVDEV
 29:            29:           
 30:           CALL DISTANCES 30:           CALL DISTANCES
 31:           DEVIATION = ABS( 100*( (NIMAGE+1)*DVEC/SEPARATION -1 ) ) 31:           DEVIATION = ABS( 100*( (NIMAGE+1)*DVEC/SEPARATION -1 ) )
 32:           AVDEV     = SUM(DEVIATION)/(NIMAGE+1) 32:           AVDEV     = SUM(DEVIATION)/(NIMAGE+1)
 33:      END SUBROUTINE IMAGEDISTRIBUTION 33:      END SUBROUTINE IMAGEDISTRIBUTION
 34:  34: 
 35:      SUBROUTINE DISTANCES     !THIS SUBROUTINE IS CALCULATING SEPARATION OF END POINTS AND 35:      SUBROUTINE DISTANCES     !THIS SUBROUTINE IS CALCULATING SEPARATION OF END POINTS AND
 36:           USE KEY,ONLY: BULKT, TWOD, BULK_BOXVEC 36:           USE KEY,ONLY: BULKT, TWOD
 37:           USE NEBDATA 37:           USE NEBDATA
 38:           USE KEYNEB,ONLY: NIMAGE 38:           USE KEYNEB,ONLY: NIMAGE
  39:           USE COMMONS,ONLY: PARAM1,PARAM2,PARAM3
 39:           USE GENRIGID ! sn402 40:           USE GENRIGID ! sn402
 40:           IMPLICIT NONE       !DISTANCES BETWEEN INDIVIDUAL IMAGES BOTH 41:           IMPLICIT NONE       !DISTANCES BETWEEN INDIVIDUAL IMAGES BOTH
 41:           DOUBLE PRECISION :: DIST_TEMP, DUMMYA(DEGFREEDOMS), DUMMYB(DEGFREEDOMS) 42:           DOUBLE PRECISION :: DIST_TEMP, DUMMYA(DEGFREEDOMS), DUMMYB(DEGFREEDOMS)
 42:           DOUBLE PRECISION :: LEFTIMAGE(DEGFREEDOMS), RIGHTIMAGE(DEGFREEDOMS) 43:           DOUBLE PRECISION :: LEFTIMAGE(DEGFREEDOMS), RIGHTIMAGE(DEGFREEDOMS)
 43:           INTEGER :: J, K 44:           INTEGER :: J, K
 44:           ! /------------------Separation(along the MEP)-------------------\ 45:           ! /------------------Separation(along the MEP)-------------------\
 45:           ! Q--------Image1------Image2---- ...----ImageNimage-------------Fin 46:           ! Q--------Image1------Image2---- ...----ImageNimage-------------Fin
 46:           ! \_DVec(1)_/ \_DVec(2)_/                       \_DVec(Nimage+1)_/ 47:           ! \_DVec(1)_/ \_DVec(2)_/                       \_DVec(Nimage+1)_/
 47:           ! sn402: added this section 48:           ! sn402: added this section
 48:           IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN 49:           IF (RIGIDINIT .AND. .NOT. ATOMRIGIDCOORDT) THEN
 50:                 LEFTIMAGE(:) = XYZ((J-1)*NOPT+1:(J-1)*NOPT+DEGFREEDOMS) 51:                 LEFTIMAGE(:) = XYZ((J-1)*NOPT+1:(J-1)*NOPT+DEGFREEDOMS)
 51:                 RIGHTIMAGE(:) = XYZ(J*NOPT+1:J*NOPT+DEGFREEDOMS) 52:                 RIGHTIMAGE(:) = XYZ(J*NOPT+1:J*NOPT+DEGFREEDOMS)
 52:                 CALL RB_DISTANCE(DIST_TEMP, LEFTIMAGE, RIGHTIMAGE, DUMMYA, DUMMYB, .FALSE.) 53:                 CALL RB_DISTANCE(DIST_TEMP, LEFTIMAGE, RIGHTIMAGE, DUMMYA, DUMMYB, .FALSE.)
 53:                 DVEC(J) = DIST_TEMP*DIST_TEMP 54:                 DVEC(J) = DIST_TEMP*DIST_TEMP
 54:             ENDDO 55:             ENDDO
 55:           ELSE IF (BULKT) THEN 56:           ELSE IF (BULKT) THEN
 56:              DO J=1,NIMAGE+1 57:              DO J=1,NIMAGE+1
 57:                 DVEC(J)=0.0D0 58:                 DVEC(J)=0.0D0
 58:                 DO K=1,NATOMS 59:                 DO K=1,NATOMS
 59:                    DVEC(J)=DVEC(J)+( XYZ(NOPT*J+3*(K-1)+1) - XYZ(NOPT*(J-1)+3*(K-1)+1) & 60:                    DVEC(J)=DVEC(J)+( XYZ(NOPT*J+3*(K-1)+1) - XYZ(NOPT*(J-1)+3*(K-1)+1) &
 60:   &                    -BULK_BOXVEC(1)*NINT((XYZ(NOPT*J+3*(K-1)+1) & 61:   &                    -PARAM1*NINT((XYZ(NOPT*J+3*(K-1)+1) - XYZ(NOPT*(J-1)+3*(K-1)+1))/PARAM1) )**2 &
 61:   &                    -XYZ(NOPT*(J-1)+3*(K-1)+1))/BULK_BOXVEC(1)) )**2 & 
 62:   &                               +( XYZ(NOPT*J+3*(K-1)+2) - XYZ(NOPT*(J-1)+3*(K-1)+2 ) & 62:   &                               +( XYZ(NOPT*J+3*(K-1)+2) - XYZ(NOPT*(J-1)+3*(K-1)+2 ) &
 63:   &                    -BULK_BOXVEC(2)*NINT((XYZ(NOPT*J+3*(K-1)+2) & 63:   &                    -PARAM2*NINT((XYZ(NOPT*J+3*(K-1)+2) - XYZ(NOPT*(J-1)+3*(K-1)+2))/PARAM2) )**2 
 64:   &                    -XYZ(NOPT*(J-1)+3*(K-1)+2))/BULK_BOXVEC(2)) )**2  
 65:                    IF (.NOT.TWOD) DVEC(J)=DVEC(J) & 64:                    IF (.NOT.TWOD) DVEC(J)=DVEC(J) &
 66:                                   +( XYZ(NOPT*J+3*(K-1)+3) - XYZ(NOPT*(J-1)+3*(K-1)+3 ) & 65:                                   +( XYZ(NOPT*J+3*(K-1)+3) - XYZ(NOPT*(J-1)+3*(K-1)+3 ) &
 67:   &                    -BULK_BOXVEC(3)*NINT((XYZ(NOPT*J+3*(K-1)+3) & 66:   &                    -PARAM3*NINT((XYZ(NOPT*J+3*(K-1)+3) - XYZ(NOPT*(J-1)+3*(K-1)+3))/PARAM3) )**2
 68:   &                    -XYZ(NOPT*(J-1)+3*(K-1)+3))/BULK_BOXVEC(3)) )**2 
 69:                 ENDDO 67:                 ENDDO
 70:              ENDDO 68:              ENDDO
 71:           ELSE 69:           ELSE
 72:              DO J=1,NIMAGE+1 70:              DO J=1,NIMAGE+1
 73:                   DVEC(J)   = SUM( ( XYZ( NOPT*J+1:NOPT*(J+1) ) - XYZ( NOPT*(J-1)+1:NOPT*J ) )**2) 71:                   DVEC(J)   = SUM( ( XYZ( NOPT*J+1:NOPT*(J+1) ) - XYZ( NOPT*(J-1)+1:NOPT*J ) )**2)
 74:              ENDDO 72:              ENDDO
 75:           ENDIF 73:           ENDIF
 76:           DVEC(1:NIMAGE+1)=SQRT(DVEC(1:NIMAGE+1)) 74:           DVEC(1:NIMAGE+1)=SQRT(DVEC(1:NIMAGE+1))
 77:           SEPARATION=SUM(DVEC(1:NIMAGE+1)) 75:           SEPARATION=SUM(DVEC(1:NIMAGE+1))
 78:  76: 
112:      END SUBROUTINE INTERNALIMAGEDISTRIBUTION110:      END SUBROUTINE INTERNALIMAGEDISTRIBUTION
113: 111: 
114: !--------------------------------------------------------------------------------------------------------112: !--------------------------------------------------------------------------------------------------------
115: 113: 
116:      SUBROUTINE MAKEIMAGE(EINITIAL,EFINAL,QQ,FINFIN) ! INTERPOLATES THE BAND114:      SUBROUTINE MAKEIMAGE(EINITIAL,EFINAL,QQ,FINFIN) ! INTERPOLATES THE BAND
117:           USE NEBDATA115:           USE NEBDATA
118:           USE SPFUNCTS116:           USE SPFUNCTS
119: !          USE KEYNEB,ONLY: NIMAGE117: !          USE KEYNEB,ONLY: NIMAGE
120:           USE KEY,ONLY: MORPHT, MSTEPS, GREATCIRCLET, GCIMAGE, GCCONV, GCUPDATE, GCSTEPS, INTEPSILON, &118:           USE KEY,ONLY: MORPHT, MSTEPS, GREATCIRCLET, GCIMAGE, GCCONV, GCUPDATE, GCSTEPS, INTEPSILON, &
121:                         BULKT, RBAAT, NTSITES, AMBERT, LOCALPERMDIST, NRBGROUP, RBGROUP, RBNINGROUP, NRBTRIES, NABT,TWOD, &119:                         BULKT, RBAAT, NTSITES, AMBERT, LOCALPERMDIST, NRBGROUP, RBGROUP, RBNINGROUP, NRBTRIES, NABT,TWOD, &
122:                         RIGIDBODY, AVOID_COLLISIONS, COLL_TOL, BULK_BOXVEC120:                         RIGIDBODY, AVOID_COLLISIONS, COLL_TOL
123:           USE INTCOMMONS, ONLY : NNZ, KD, NINTC, DESMINT, DIHINFO, PREVDIH, BACKTCUTOFF, INTERPBACKTCUT, MINBACKTCUT, &121:           USE INTCOMMONS, ONLY : NNZ, KD, NINTC, DESMINT, DIHINFO, PREVDIH, BACKTCUTOFF, INTERPBACKTCUT, MINBACKTCUT, &
124:                                  CHICDNEB122:                                  CHICDNEB
125:           USE COMMONS, ONLY : ZSYM, NRBSITES, DEBUG123:           USE COMMONS, ONLY : ZSYM, PARAM1,PARAM2,PARAM3, NRBSITES, DEBUG
126:           USE MODCHARMM, ONLY : CHRMMT, ICINTERPT124:           USE MODCHARMM, ONLY : CHRMMT, ICINTERPT
127:           USE MODAMBER9, ONLY: AMBERICT, AMBICDNEBT, NICTOT125:           USE MODAMBER9, ONLY: AMBERICT, AMBICDNEBT, NICTOT
128:           USE PORFUNCS 126:           USE PORFUNCS 
129:           USE KEYNEB,ONLY: NIMAGE,XYZFILE,RBXYZFILE127:           USE KEYNEB,ONLY: NIMAGE,XYZFILE,RBXYZFILE
130:           USE KEY, ONLY: FILTH,FILTHSTR, MULTISITEPYT, ALIGNRBST, N_TO_ALIGN, TO_ALIGN !sn402128:           USE KEY, ONLY: FILTH,FILTHSTR, MULTISITEPYT, ALIGNRBST, N_TO_ALIGN, TO_ALIGN !sn402
131:           USE GENRIGID, ONLY: RIGIDMOLECULEST, DEGFREEDOMS, NRIGIDBODY, GENRIGID_IMAGE_CTORIGID, GENRIGID_IMAGE_RIGIDTOC !sn402129:           USE GENRIGID, ONLY: RIGIDMOLECULEST, DEGFREEDOMS, NRIGIDBODY, GENRIGID_IMAGE_CTORIGID, GENRIGID_IMAGE_RIGIDTOC !sn402
132:           IMPLICIT NONE130:           IMPLICIT NONE
133:           INTEGER :: I, J1, ITDONE, INTERVAL, NDONE, J2, J, K, ISTAT, J3, J4, J5, NDUMMY, J6, J7, JWORST2, JWORST3131:           INTEGER :: I, J1, ITDONE, INTERVAL, NDONE, J2, J, K, ISTAT, J3, J4, J5, NDUMMY, J6, J7, JWORST2, JWORST3
134:           DOUBLE PRECISION,ALLOCATABLE :: DELTAX(:)132:           DOUBLE PRECISION,ALLOCATABLE :: DELTAX(:)
135: !          DOUBLE PRECISION, ALLOCATABLE :: QTN(:,:), PTN(:,:)133: !          DOUBLE PRECISION, ALLOCATABLE :: QTN(:,:), PTN(:,:)
259:              ALLOCATE(DELTAX(NOPT))257:              ALLOCATE(DELTAX(NOPT))
260: 258: 
261:             ! DELTAX is the change in each coordinate between each consecutive pair of images259:             ! DELTAX is the change in each coordinate between each consecutive pair of images
262:             ! Note: we are performing some unnecessary calculations here because we260:             ! Note: we are performing some unnecessary calculations here because we
263:             ! compute DELTAX for the angle-axis coordinates as well as the CoM coordinates, but we never261:             ! compute DELTAX for the angle-axis coordinates as well as the CoM coordinates, but we never
264:             ! use the former. I hope this won't be a big problem for performance so I've left it in to keep262:             ! use the former. I hope this won't be a big problem for performance so I've left it in to keep
265:             ! the code tidy.263:             ! the code tidy.
266:             IF (BULKT) THEN264:             IF (BULKT) THEN
267:                 DO K=1,NATOMS265:                 DO K=1,NATOMS
268:                    DELTAX(3*(K-1)+1)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1) &266:                    DELTAX(3*(K-1)+1)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1) &
269:   &                    -BULK_BOXVEC(1)*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1))/BULK_BOXVEC(1))267:   &                    -PARAM1*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1))/PARAM1)
270:                    DELTAX(3*(K-1)+2)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+2) - XYZ(3*(K-1)+2) &268:                    DELTAX(3*(K-1)+2)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+2) - XYZ(3*(K-1)+2) &
271:   &                    -BULK_BOXVEC(2)*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+2) - XYZ(3*(K-1)+2))/BULK_BOXVEC(2))269:   &                    -PARAM2*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+2) - XYZ(3*(K-1)+2))/PARAM2)
272:                    IF (.NOT.TWOD) DELTAX(3*(K-1)+3)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+3) - XYZ(3*(K-1)+3) &270:                    IF (.NOT.TWOD) DELTAX(3*(K-1)+3)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+3) - XYZ(3*(K-1)+3) &
273:   &                    -BULK_BOXVEC(3)*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+3) - XYZ(3*(K-1)+3))/BULK_BOXVEC(3))271:   &                    -PARAM3*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+3) - XYZ(3*(K-1)+3))/PARAM3)
274:                 ENDDO272:                 ENDDO
275:                 DELTAX(1:NOPT)=DELTAX(1:NOPT)/(NIMAGE+1)273:                 DELTAX(1:NOPT)=DELTAX(1:NOPT)/(NIMAGE+1)
276:             ELSE274:             ELSE
277:                 DELTAX(1:NOPT) = ( XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2) ) - XYZ(1:NOPT) )/(NIMAGE+1)275:                 DELTAX(1:NOPT) = ( XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2) ) - XYZ(1:NOPT) )/(NIMAGE+1)
278:             ENDIF276:             ENDIF
279:             277:             
280:             ! Linear interpolation of the centre of mass coordinates for all images278:             ! Linear interpolation of the centre of mass coordinates for all images
281:              DO I = 2, NIMAGE+1279:              DO I = 2, NIMAGE+1
282:                 J = NOPT*(I-1)280:                 J = NOPT*(I-1)
283:                 IF(RIGIDMOLECULEST) THEN281:                 IF(RIGIDMOLECULEST) THEN
373: 371: 
374:             DEALLOCATE(DELTAX)372:             DEALLOCATE(DELTAX)
375: 373: 
376:           ELSE   ! End of "else if(RBAAT .OR. RIGIDMOLECULEST)"374:           ELSE   ! End of "else if(RBAAT .OR. RIGIDMOLECULEST)"
377: 375: 
378:              ALLOCATE(DELTAX(NOPT))376:              ALLOCATE(DELTAX(NOPT))
379: 377: 
380:              IF (BULKT) THEN378:              IF (BULKT) THEN
381:                 DO K=1,NATOMS379:                 DO K=1,NATOMS
382:                    DELTAX(3*(K-1)+1)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1) &380:                    DELTAX(3*(K-1)+1)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1) &
383:   &                    -BULK_BOXVEC(1)*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1))/BULK_BOXVEC(1)) 381:   &                    -PARAM1*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+1) - XYZ(3*(K-1)+1))/PARAM1) 
384:                    DELTAX(3*(K-1)+2)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+2) - XYZ(3*(K-1)+2) &382:                    DELTAX(3*(K-1)+2)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+2) - XYZ(3*(K-1)+2) &
385:   &                    -BULK_BOXVEC(2)*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+2) - XYZ(3*(K-1)+2))/BULK_BOXVEC(2))383:   &                    -PARAM2*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+2) - XYZ(3*(K-1)+2))/PARAM2)
386:                    IF (.NOT.TWOD) DELTAX(3*(K-1)+3)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+3) - XYZ(3*(K-1)+3) &384:                    IF (.NOT.TWOD) DELTAX(3*(K-1)+3)=XYZ(NOPT*(NIMAGE+1)+3*(K-1)+3) - XYZ(3*(K-1)+3) &
387:   &                    -BULK_BOXVEC(3)*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+3) - XYZ(3*(K-1)+3))/BULK_BOXVEC(3))385:   &                    -PARAM3*NINT((XYZ(NOPT*(NIMAGE+1)+3*(K-1)+3) - XYZ(3*(K-1)+3))/PARAM3)
388:                 ENDDO386:                 ENDDO
389:                 DELTAX(1:NOPT)=DELTAX(1:NOPT)/(NIMAGE+1)387:                 DELTAX(1:NOPT)=DELTAX(1:NOPT)/(NIMAGE+1)
390:              ELSE388:              ELSE
391:                 DELTAX(1:NOPT) = ( XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2) ) - XYZ(1:NOPT) )/(NIMAGE+1)389:                 DELTAX(1:NOPT) = ( XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2) ) - XYZ(1:NOPT) )/(NIMAGE+1)
392:              ENDIF390:              ENDIF
393:    391:    
394: ! bs360: Interpolation using dihedrals.392: ! bs360: Interpolation using dihedrals.
395:              IF (CHRMMT.AND.CHICDNEB.AND.ICINTERPT) THEN393:              IF (CHRMMT.AND.CHICDNEB.AND.ICINTERPT) THEN
396:                 QS(1:NOPT)=XYZ(1:NOPT)         394:                 QS(1:NOPT)=XYZ(1:NOPT)         
397:                 QF(1:NOPT)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2))395:                 QF(1:NOPT)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+2))
498:            QF(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)496:            QF(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)
499: !497: !
500: !  Now change the whole XYZ array to be based on the permutational isomer that498: !  Now change the whole XYZ array to be based on the permutational isomer that
501: !  minimises the overall distance. The difference is in the initial straight line499: !  minimises the overall distance. The difference is in the initial straight line
502: !  guess, which may be for different permutational isomers.500: !  guess, which may be for different permutational isomers.
503: !  Use the permutation that gives the lowest value for the highest image.501: !  Use the permutation that gives the lowest value for the highest image.
504: !502: !
505:            QS2(1:3*NATOMS)=XYZ(1:NOPT)503:            QS2(1:3*NATOMS)=XYZ(1:NOPT)
506:            QF2(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)504:            QF2(1:3*NATOMS)=XYZ(NOPT*(NIMAGE+1)+1:NOPT*(NIMAGE+1)+NOPT)
507:            CALL MINPERMDIST(QS2,QF2,NATOMS, &505:            CALL MINPERMDIST(QS2,QF2,NATOMS, &
508:   &                         DEBUG,BULK_BOXVEC(1),BULK_BOXVEC(2),BULK_BOXVEC(3),BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)506:   &                         DEBUG,PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
509:            EWORST2=-HUGE(1.0D0)507:            EWORST2=-HUGE(1.0D0)
510:            DO J1=1,NIMAGE508:            DO J1=1,NIMAGE
511:               COORDS(1:3*NATOMS)=(J1*QF2(1:3*NATOMS)+(NIMAGE+1-J1)*QS2(1:3*NATOMS))/(NIMAGE+1)509:               COORDS(1:3*NATOMS)=(J1*QF2(1:3*NATOMS)+(NIMAGE+1-J1)*QS2(1:3*NATOMS))/(NIMAGE+1)
512:               CALL POTENTIAL(COORDS,E1,LGDUMMY,.FALSE.,.FALSE.,RMSDUMMY,.FALSE.,.FALSE.)510:               CALL POTENTIAL(COORDS,E1,LGDUMMY,.FALSE.,.FALSE.,RMSDUMMY,.FALSE.,.FALSE.)
513:               LEIMAGE2(J1)=E1511:               LEIMAGE2(J1)=E1
514:               IF (E1.GT.EWORST2) THEN512:               IF (E1.GT.EWORST2) THEN
515:                  JWORST2=J1513:                  JWORST2=J1
516:                  EWORST2=E1514:                  EWORST2=E1
517:               ENDIF515:               ENDIF
518: !             IF (DEBUG) PRINT '(A,I6,A,G20.10)',' nnutils> alternative interpolation, image ',J1,' energy=',E1516: !             IF (DEBUG) PRINT '(A,I6,A,G20.10)',' nnutils> alternative interpolation, image ',J1,' energy=',E1
1058:     DOUBLE PRECISION, INTENT(OUT) :: EWORST1056:     DOUBLE PRECISION, INTENT(OUT) :: EWORST
1059:     INTEGER, INTENT(OUT) :: JWORST ! Energy and index of the highest-energy index, respectively1057:     INTEGER, INTENT(OUT) :: JWORST ! Energy and index of the highest-energy index, respectively
1060:     INTEGER :: J11058:     INTEGER :: J1
1061:     DOUBLE PRECISION :: E1, RMSDUMMY, LGDUMMY(NOPT), COORDS(NOPT)1059:     DOUBLE PRECISION :: E1, RMSDUMMY, LGDUMMY(NOPT), COORDS(NOPT)
1062: 1060: 
1063:     EWORST=-HUGE(1.0D0)1061:     EWORST=-HUGE(1.0D0)
1064:     1062:     
1065:     DO J1=1,NIMAGE1063:     DO J1=1,NIMAGE
1066: 1064: 
1067:         COORDS = BAND(NOPT*J1+1:NOPT*(J1+1))1065:         COORDS = BAND(NOPT*J1+1:NOPT*(J1+1))
 1066:         write(*,*) "Image coords"
 1067:         write(*,*) COORDS
1068:         ! Determine the energy of this image1068:         ! Determine the energy of this image
1069:         CALL POTENTIAL(COORDS,E1,LGDUMMY,.FALSE.,.FALSE.,RMSDUMMY,.FALSE.,.FALSE.)1069:         CALL POTENTIAL(COORDS,E1,LGDUMMY,.FALSE.,.FALSE.,RMSDUMMY,.FALSE.,.FALSE.)
1070: 1070: 
1071:         IF (E1.GT.EWORST) THEN1071:         IF (E1.GT.EWORST) THEN
1072:             JWORST=J11072:             JWORST=J1
1073:             EWORST=E11073:             EWORST=E1
1074:         ENDIF1074:         ENDIF
1075:         IF (DEBUG) PRINT '(A,I6,A,G20.10)',' nnutils> image ',J1,' energy=',E11075:         IF (DEBUG) PRINT '(A,I6,A,G20.10)',' nnutils> image ',J1,' energy=',E1
1076:     ENDDO1076:     ENDDO
1077:     IF (DEBUG) PRINT '(A,I6,A,G20.10)',' nnutils> highest image is ',JWORST,' with energy=',EWORST 1077:     IF (DEBUG) PRINT '(A,I6,A,G20.10)',' nnutils> highest image is ',JWORST,' with energy=',EWORST 
1187: 1187: 
1188: !         IF (.NOT.PRESENT(UNITIN)) THEN1188: !         IF (.NOT.PRESENT(UNITIN)) THEN
1189:                CLOSE(UNIT)1189:                CLOSE(UNIT)
1190: !         ENDIF1190: !         ENDIF
1191:           PRINT *, 'writeprofile> NEB profile was saved to file "'//trim(filename)//'"'1191:           PRINT *, 'writeprofile> NEB profile was saved to file "'//trim(filename)//'"'
1192:      END SUBROUTINE WRITEPROFILE1192:      END SUBROUTINE WRITEPROFILE
1193: 1193: 
1194:      SUBROUTINE RWG(WHAT,GUESS,NITER)1194:      SUBROUTINE RWG(WHAT,GUESS,NITER)
1195:           USE PORFUNCS1195:           USE PORFUNCS
1196:           USE KEY,ONLY: FILTH,FILTHSTR,UNRST,STOCKT,AMHT,SEQ,NUMGLY,STOCKAAT, RBAAT,NTSITES, GTHOMSONT, NGTHORI, PERMGUESS, &1196:           USE KEY,ONLY: FILTH,FILTHSTR,UNRST,STOCKT,AMHT,SEQ,NUMGLY,STOCKAAT, RBAAT,NTSITES, GTHOMSONT, NGTHORI, PERMGUESS, &
1197:   &                     BULKT, BULK_BOXVEC, NRBTRIES, NABT,TWOD, RIGIDBODY, VARIABLES1197:   &                     BULKT, NRBTRIES, NABT,TWOD, RIGIDBODY, VARIABLES
1198:           USE COMMONS, ONLY: ZSYM, NRBSITES, NRBSITES, DEBUG1198:           USE COMMONS, ONLY: ZSYM, NRBSITES, PARAM1,PARAM2,PARAM3, NRBSITES, DEBUG
1199:           USE INTCOMMONS, ONLY : DESMINT1199:           USE INTCOMMONS, ONLY : DESMINT
1200:           USE NEBDATA1200:           USE NEBDATA
1201:           USE AMHGLOBALS, ONLY : NMRES1201:           USE AMHGLOBALS, ONLY : NMRES
1202:           USE KEYNEB,ONLY: NIMAGE,XYZFILE,RBXYZFILE,GUESSFILE1202:           USE KEYNEB,ONLY: NIMAGE,XYZFILE,RBXYZFILE,GUESSFILE
1203:           USE GENRIGID1203:           USE GENRIGID
1204:           IMPLICIT NONE1204:           IMPLICIT NONE
1205: 1205: 
1206:           CHARACTER,INTENT(IN) :: WHAT1206:           CHARACTER,INTENT(IN) :: WHAT
1207:           LOGICAL,INTENT(IN) :: GUESS1207:           LOGICAL,INTENT(IN) :: GUESS
1208:           INTEGER,INTENT(IN) :: NITER1208:           INTEGER,INTENT(IN) :: NITER
1356: !              ENDIF1356: !              ENDIF
1357:                DO J2=1,NIMAGE+21357:                DO J2=1,NIMAGE+2
1358: !1358: !
1359: ! Here we skip two lines, allowing for the second line to be blank.1359: ! Here we skip two lines, allowing for the second line to be blank.
1360: !1360: !
1361:                   READ(993,*)1361:                   READ(993,*)
1362:                   READ(993,*)1362:                   READ(993,*)
1363:                   READ(993,*) (DUMMYS,XYZ( (J2-1)*NOPT+J1),XYZ((J2-1)*NOPT+J1+1),XYZ((J2-1)*NOPT+J1+2),J1=1,NOPT,3)1363:                   READ(993,*) (DUMMYS,XYZ( (J2-1)*NOPT+J1),XYZ((J2-1)*NOPT+J1+1),XYZ((J2-1)*NOPT+J1+2),J1=1,NOPT,3)
1364:                   IF (PERMGUESS.AND.(J2.GT.1)) THEN 1364:                   IF (PERMGUESS.AND.(J2.GT.1)) THEN 
1365:                      CALL MINPERMDIST(XYZ((J2-2)*NOPT+1:(J2-1)*NOPT),XYZ((J2-1)*NOPT+1:J2*NOPT),NATOMS,DEBUG, &1365:                      CALL MINPERMDIST(XYZ((J2-2)*NOPT+1:(J2-1)*NOPT),XYZ((J2-1)*NOPT+1:J2*NOPT),NATOMS,DEBUG, &
1366:   &                    BULK_BOXVEC(1),BULK_BOXVEC(2),BULK_BOXVEC(3),BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)1366:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
1367:                   ENDIF1367:                   ENDIF
1368:                ENDDO1368:                ENDDO
1369:                PRINT '(A,I6,A)','nnutils> ',NIMAGE+2,' images read from file ' // TRIM(ADJUSTL(GUESSFILE))1369:                PRINT '(A,I6,A)','nnutils> ',NIMAGE+2,' images read from file ' // TRIM(ADJUSTL(GUESSFILE))
1370:           END SELECT1370:           END SELECT
1371:           CLOSE(UNIT=993)1371:           CLOSE(UNIT=993)
1372:      END SUBROUTINE RWG1372:      END SUBROUTINE RWG
1373: 1373: 
1374:      SUBROUTINE SAVEBANDCOORD1374:      SUBROUTINE SAVEBANDCOORD
1375:           USE NEBDATA,ONLY: XYZ,NOPT1375:           USE NEBDATA,ONLY: XYZ,NOPT
1376:           USE KEYNEB,ONLY:NIMAGE,PTSFILE1376:           USE KEYNEB,ONLY:NIMAGE,PTSFILE


r30520/potential.f 2016-05-26 15:30:15.212139054 +0100 r30519/potential.f 2016-05-26 15:30:17.444170191 +0100
3509:             CALL USERPOT_POTENTIALHESS(3*NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST,HESS)3509:             CALL USERPOT_POTENTIALHESS(3*NATOMS,COORDS,VNEW,ENERGY,GTEST,SSTEST,HESS)
3510:             IF (PTEST) THEN3510:             IF (PTEST) THEN
3511:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' meV'3511:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' meV'
3512:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' meV'3512:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' meV'
3513:             ENDIF3513:             ENDIF
3514: 3514: 
3515:             ELSE IF (STEALTHYT) THEN3515:             ELSE IF (STEALTHYT) THEN
3516:                CALL STEALTHY(COORDS, VNEW, ENERGY, GTEST, SSTEST)3516:                CALL STEALTHY(COORDS, VNEW, ENERGY, GTEST, SSTEST)
3517:  3517:  
3518:                !DIFF=1.0D-63518:                !DIFF=1.0D-6
3519:                !OPEN(UNIT=1500, FILE="grad-test")3519:                !PRINT*,'ANALYTIC AND NUMERICAL GRADIENTS:NATOMS=',NATOMS
3520:                !WRITE(1500,*) 'ANALYTIC AND NUMERICAL GRADIENTS:NATOMS=',NATOMS 
3521:                !DO J1=1, 3*NATOMS3520:                !DO J1=1, 3*NATOMS
3522:                !   COORDS(J1)=COORDS(J1)+DIFF3521:                !   COORDS(J1)=COORDS(J1)+DIFF
3523:                !   CALL STEALTHY(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)3522:                !   CALL STEALTHY(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)
3524:                !   COORDS(J1)=COORDS(J1)-2.0D0*DIFF3523:                !   COORDS(J1)=COORDS(J1)-2.0D0*DIFF
3525:                !   CALL STEALTHY(COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.)3524:                !   CALL STEALTHY(COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.)
3526:                !   COORDS(J1)=COORDS(J1)+DIFF               3525:                !   COORDS(J1)=COORDS(J1)+DIFF
3527:                !   IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*ABS((VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.1.0D-3)) THEN3526:                !
3528:                !      WRITE(1500,'(I5,2G20.10,A)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF), '  X'3527:                !  IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*ABS((VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.1.0D-3)) THEN
 3528:                !      WRITE(*,'(I5,2G20.10,A)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF), '  X'
3529:                !   ELSE3529:                !   ELSE
3530:                !      WRITE(1500,'(I5,2G20.10,A)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)3530:                !      WRITE(*,'(I5,2G20.10,A)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)
3531:                !   END IF3531:                !   END IF
3532:                !END DO3532:                !END DO
3533:   3533:   
3534:                !CALL STEALTHY(COORDS, VNEW, ENERGY, .TRUE., .TRUE.)3534:                !CALL STEALTHY(COORDS, VNEW, ENERGY, .TRUE., .TRUE.)
3535:                !OPEN(UNIT=1499, FILE="hessian-matrix")3535:                !OPEN(UNIT=1499, FILE="hessian-matrix")
3536:                !WRITE(1499,*) 'ANALYTIC AND NUMERICAL SECOND DERIVATIVES:'3536:                !WRITE(1499,*),'ANALYTIC AND NUMERICAL SECOND DERIVATIVES:'
3537:                !DO J1=1, 3*NATOMS3537:                !DO J1=1, 3*NATOMS
3538:                !   COORDS(J1)=COORDS(J1)+DIFF3538:                !   COORDS(J1)=COORDS(J1)+DIFF
3539:                !   CALL STEALTHY(COORDS,VPLUS,EPLUS,.TRUE.,.FALSE.)3539:                !   CALL STEALTHY(COORDS,VPLUS,EPLUS,.TRUE.,.FALSE.)
3540:                !   COORDS(J1)=COORDS(J1)-2.0D0*DIFF3540:                !   COORDS(J1)=COORDS(J1)-2.0D0*DIFF
3541:                !   CALL STEALTHY(COORDS,VMINUS,EMINUS,.TRUE.,.FALSE.)3541:                !   CALL STEALTHY(COORDS,VMINUS,EMINUS,.TRUE.,.FALSE.)
3542:                !   COORDS(J1)=COORDS(J1)+DIFF3542:                !   COORDS(J1)=COORDS(J1)+DIFF
3543:                !   DO J2=1, 3*NATOMS3543:                !   DO J2=1, 3*NATOMS
3544:                !      IF (ABS(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)).GT.(1.0D-5*HESS(J1,J2))) THEN3544:                !      IF (ABS(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)).GT.(1.0D-5*HESS(J1,J2))) THEN
3545:                !         WRITE(1499,'(2I5,2F20.15,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'  X'3545:                !         WRITE(1499,'(2I5,2F20.15,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'  X'
3546:                !      ELSE3546:                !      ELSE


r30520/stealthy.f90 2016-05-26 15:30:15.436142178 +0100 r30519/stealthy.f90 2016-05-26 15:30:17.676173427 +0100
  5:     IMPLICIT NONE  5:     IMPLICIT NONE
  6:     LOGICAL GTEST,STEST  6:     LOGICAL GTEST,STEST
  7:     INTEGER M, M1, M2, M3, J1, J2, I1, I2, L1, L2, KI  7:     INTEGER M, M1, M2, M3, J1, J2, I1, I2, L1, L2, KI
  8:     DOUBLE PRECISION XS(3*NATOMS), KX(3), KY(3), KZ(3), KN(3), TIDU(3*NATOMS), &  8:     DOUBLE PRECISION XS(3*NATOMS), KX(3), KY(3), KZ(3), KN(3), TIDU(3*NATOMS), &
  9:                      CKR(NATOMS), SKR(NATOMS)  9:                      CKR(NATOMS), SKR(NATOMS)
 10:     DOUBLE PRECISION PI, KR ,C, S, ENERGY, VF, IM, KM 10:     DOUBLE PRECISION PI, KR ,C, S, ENERGY, VF, IM, KM
 11:  11: 
 12:     IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS, 3*NATOMS)) 12:     IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS, 3*NATOMS))
 13:  13: 
 14: !   initialize 14: !   initialize
 15:     PI=3.1415926535897932384626 15:     PI=3.141592653589793
 16:     KM=KLIM*KLIM 16:     KM=KLIM*KLIM
 17:     KI=INT(KLIM) 17:     KI=DINT(KLIM)
 18:     ENERGY=0 18:     ENERGY=0
 19:     STM=0 
 20:  19: 
 21:     IF (GTEST) THEN 20:     IF (GTEST) THEN
 22:        TIDU=0 21:        TIDU=0
 23:     END IF 22:     END IF
 24:  23: 
 25:     IF (STEST) THEN 24:     IF (STEST) THEN
 26:        HESS=0 25:        HESS=0
 27:     END IF 26:     END IF
 28:  27: 
 29: !   get the initial box 28: !   get the initial box
 30:  29: 
 31: !   calculate the wave vector 30: !   calculate the wave vector
 32:     KX(1)=360/BULK_BOXVEC(1) 31:     KX(1)=2*PI/BULK_BOXVEC(1)
 33:     KX(2)=0 32:     KX(2)=0
 34:     KX(3)=0 33:     KX(3)=0
 35:     KY(1)=0 34:     KY(1)=0
 36:     KY(2)=360/BULK_BOXVEC(2) 35:     KY(2)=2*PI/BULK_BOXVEC(2)
 37:     KY(3)=0 36:     KY(3)=0
 38:     KZ(1)=0 37:     KZ(1)=0
 39:     KZ(2)=0 38:     KZ(2)=0
 40:     KZ(3)=360/BULK_BOXVEC(3) 39:     KZ(3)=2*PI/BULK_BOXVEC(3)
 41:     VF=BULK_BOXVEC(1)*BULK_BOXVEC(2)*BULK_BOXVEC(3) 40:     VF=BULK_BOXVEC(1)*BULK_BOXVEC(2)*BULK_BOXVEC(3)
 42: !   WRITE(*,*) KI(1), KI(2), KI(3) 41: !   WRITE(*,*) KI(1), KI(2), KI(3)
 43:  42: 
 44: !   periodic boundary condition 43: !   periodic boundary condition
 45:     IF (BULKT.AND. .NOT.REBOXT) THEN 44:     IF (BULKT) THEN
 46:        DO J1=1, NATOMS 45:        DO J1=1, NATOMS
 47:           J2=3*(J1-1) 46:           J2=3*(J1-1)
 48:           XS(J2+1)=XS(J2+1) - BULK_BOXVEC(1)*DNINT(XS(J2+1)/BULK_BOXVEC(1)) 47:           XS(J2+1)=XS(J2+1) - BULK_BOXVEC(1)*DNINT(XS(J2+1)/BULK_BOXVEC(1))
 49:           XS(J2+2)=XS(J2+2) - BULK_BOXVEC(2)*DNINT(XS(J2+2)/BULK_BOXVEC(2)) 48:           XS(J2+2)=XS(J2+2) - BULK_BOXVEC(2)*DNINT(XS(J2+2)/BULK_BOXVEC(2))
 50:           XS(J2+3)=XS(J2+3) - BULK_BOXVEC(3)*DNINT(XS(J2+3)/BULK_BOXVEC(3)) 49:           XS(J2+3)=XS(J2+3) - BULK_BOXVEC(3)*DNINT(XS(J2+3)/BULK_BOXVEC(3))
 51:        END DO 50:        END DO
 52:     END IF 51:     END IF
 53: !    write(*,*) XS 
 54:   52:  
 55: !   sum of k 53: !   sum of k
 56:     DO M1=0, KI 54:     DO M1=0, KI
 57:        DO M2=-KI, KI 55:        DO M2=-KI, KI
 58:           DO M3=-KI, KI 56:           DO M3=-KI, KI
 59:  57: 
 60:              IF (M1==0.AND.M2>=0.AND.M3<0) CYCLE 58:              IF (M1==0.AND.M2>=0.AND.M3<0) CYCLE
 61:              IF (M1==0.AND.M2<=0.AND.M3<=0) CYCLE 59:              IF (M1==0.AND.M2<=0.AND.M3<=0) CYCLE
 62:  60: 
 63:              M=M1*M1+M2*M2+M3*M3 61:              M=M1*M1+M2*M2+M3*M3
 66:              C=0 64:              C=0
 67:              S=0 65:              S=0
 68:              KN(1)=KX(1)*M1 + KY(1)*M2 + KZ(1)*M3 66:              KN(1)=KX(1)*M1 + KY(1)*M2 + KZ(1)*M3
 69:              KN(2)=KX(2)*M1 + KY(2)*M2 + KZ(2)*M3 67:              KN(2)=KX(2)*M1 + KY(2)*M2 + KZ(2)*M3
 70:              KN(3)=KX(3)*M1 + KY(3)*M2 + KZ(3)*M3 68:              KN(3)=KX(3)*M1 + KY(3)*M2 + KZ(3)*M3
 71:  69: 
 72: !            calculation of Re[n(k)] and Im[n(k)] 70: !            calculation of Re[n(k)] and Im[n(k)]
 73:              DO I1=1, NATOMS 71:              DO I1=1, NATOMS
 74:                 I2=3*(I1-1) 72:                 I2=3*(I1-1)
 75:                 KR=KN(1)*XS(I2+1) + KN(2)*XS(I2+2) + KN(3)*XS(I2+3) 73:                 KR=KN(1)*XS(I2+1) + KN(2)*XS(I2+2) + KN(3)*XS(I2+3)
 76:                 CKR(I1)=COSD(KR) 74:                 CKR(I1)=COS(KR)
 77:                 SKR(I1)=SIND(KR) 75:                 SKR(I1)=SIN(KR)
 78:                 !if (M1==0.AND.M2==-3) write(*,'(T1,4(I,1X),F,1X,F)') M1, M2, M3, I1, CKR(I1), SKR(I1) 
 79:  
 80:                 C=C+CKR(I1) 76:                 C=C+CKR(I1)
 81:                 S=S+SKR(I1) 77:                 S=S+SKR(I1)
 82: !                if (M1==0.AND.M2==-3) write(*,'(T1,F,1X,F)') C, S 
 83:              END DO 78:              END DO
 84:              !write(*,*) M1, M2, M3 
 85:              !write(*,*) CKR(1:6) 
 86:              !write(*,*) SKR(1:6) 
 87:              !if (M1==0.AND.M2==-8) then 
 88:              !   write(*,*) CKR(1:6) 
 89:              !   write(*,*) SKR(1:6) 
 90:              !end if 
 91: !             write(*,'(T1,3(I,1X),F,1X,F)') M1, M2, M3, C, S 
 92:  79: 
 93: !            E=n(k)^2/Vf 80: !            E=n(k)^2/Vf
 94:              ENERGY=ENERGY+ SCA*(C*C+S*S)/VF 81:              ENERGY=ENERGY+ SCA*(C*C+S*S)/VF
 95:  82: 
 96: !            gradient 83: !            gradient
 97:              IF (GTEST) THEN 84:              IF (GTEST) THEN
 98:                 DO L1=1, NATOMS 85:                 DO L1=1, NATOMS
 99:                    L2=3*(L1-1) 86:                    L2=3*(L1-1)
100:                    IM=( C*SKR(L1) - S*CKR(L1) )*SCA 87:                    IM=( C*SKR(L1) - S*CKR(L1) )*SCA
101:                    TIDU(L2+1)=TIDU(L2+1) - 2*KN(1)*IM/VF 88:                    TIDU(L2+1)=TIDU(L2+1) - 2*KN(1)*IM/VF
102:                    TIDU(L2+2)=TIDU(L2+2) - 2*KN(2)*IM/VF 89:                    TIDU(L2+2)=TIDU(L2+2) - 2*KN(2)*IM/VF
103:                    TIDU(L2+3)=TIDU(L2+3) - 2*KN(3)*IM/VF 90:                    TIDU(L2+3)=TIDU(L2+3) - 2*KN(3)*IM/VF
104:                 END DO 91:                 END DO
105:              END IF 92:              END IF
106:  93: 
107: !            hessian matrix 94: !            hessian matrix
108:              IF (STEST) THEN 95:              IF (STEST) THEN
109:                 CALL STHESS(XS, KN, VF, CKR, C, SKR, S) 96:                 CALL STHESS(XS, KN, VF, CKR, C, SKR, S)
110:              END IF 97:              END IF
111:  98: 
112:              STM=STM+1 
113:  
114:           END DO 99:           END DO
115:        END DO100:        END DO
116:     END DO101:     END DO
117: 102:    
118:     IF (STEST.AND.STEALTV) THEN 
119:        CALL STEALTHY_EIGEN_TEST(ENERGY, TIDU) 
120:     ENDIF 
121:  
122: !   change the box back to cubic103: !   change the box back to cubic
123: 104: 
124: END SUBROUTINE105: END SUBROUTINE
125: 106: 
126: SUBROUTINE STHESS(X, K, V, CK, CSUM, SK, SSUM)107: SUBROUTINE STHESS(X, K, V, CK, CSUM, SK, SSUM)
127:    USE COMMONS108:    USE COMMONS
128:    USE KEY109:    USE KEY
129:    USE MODHESS110:    USE MODHESS
130:    IMPLICIT NONE111:    IMPLICIT NONE
131:    INTEGER I1, I2, J1, J2112:    INTEGER I1, I2, J1, J2
147:       HESS(I2+3, I2+3)=HESS(I2+3, I2+3) - 2*K(3)*K(3)/V * CF128:       HESS(I2+3, I2+3)=HESS(I2+3, I2+3) - 2*K(3)*K(3)/V * CF
148: 129: 
149:       DO J1=I1+1, NATOMS130:       DO J1=I1+1, NATOMS
150:          J2=3*(J1-1)131:          J2=3*(J1-1)
151: 132: 
152:          RIJ(1)=X(I2+1) - X(J2+1)133:          RIJ(1)=X(I2+1) - X(J2+1)
153:          RIJ(2)=X(I2+2) - X(J2+2)134:          RIJ(2)=X(I2+2) - X(J2+2)
154:          RIJ(3)=X(I2+3) - X(J2+3)135:          RIJ(3)=X(I2+3) - X(J2+3)
155: 136: 
156:          KRIJ=K(1)*RIJ(1) + K(2)*RIJ(2) + K(3)*RIJ(3)137:          KRIJ=K(1)*RIJ(1) + K(2)*RIJ(2) + K(3)*RIJ(3)
157:          CSKR= COSD(KRIJ)*SCA138:          CSKR= COS(KRIJ)*SCA
158: 139: 
159:          HESS(I2+1, J2+1)=HESS(I2+1, J2+1) + 2*K(1)*K(1)/V * CSKR140:          HESS(I2+1, J2+1)=HESS(I2+1, J2+1) + 2*K(1)*K(1)/V * CSKR
160:          HESS(I2+1, J2+2)=HESS(I2+1, J2+2) + 2*K(1)*K(2)/V * CSKR141:          HESS(I2+1, J2+2)=HESS(I2+1, J2+2) + 2*K(1)*K(2)/V * CSKR
161:          HESS(I2+1, J2+3)=HESS(I2+1, J2+3) + 2*K(1)*K(3)/V * CSKR142:          HESS(I2+1, J2+3)=HESS(I2+1, J2+3) + 2*K(1)*K(3)/V * CSKR
162: 143: 
163:          HESS(I2+2, J2+1)=HESS(I2+2, J2+1) + 2*K(2)*K(1)/V * CSKR144:          HESS(I2+2, J2+1)=HESS(I2+2, J2+1) + 2*K(2)*K(1)/V * CSKR
164:          HESS(I2+2, J2+2)=HESS(I2+2, J2+2) + 2*K(2)*K(2)/V * CSKR145:          HESS(I2+2, J2+2)=HESS(I2+2, J2+2) + 2*K(2)*K(2)/V * CSKR
165:          HESS(I2+2, J2+3)=HESS(I2+2, J2+3) + 2*K(2)*K(3)/V * CSKR146:          HESS(I2+2, J2+3)=HESS(I2+2, J2+3) + 2*K(2)*K(3)/V * CSKR
166: 147: 
167:          HESS(I2+3, J2+1)=HESS(I2+3, J2+1) + 2*K(3)*K(1)/V * CSKR148:          HESS(I2+3, J2+1)=HESS(I2+3, J2+1) + 2*K(3)*K(1)/V * CSKR
190:          HESS(J2+3, I2+2)=HESS(I2+2, J2+3)171:          HESS(J2+3, I2+2)=HESS(I2+2, J2+3)
191: 172: 
192:          HESS(J2+1, I2+3)=HESS(I2+3, J2+1)173:          HESS(J2+1, I2+3)=HESS(I2+3, J2+1)
193:          HESS(J2+2, I2+3)=HESS(I2+3, J2+2)174:          HESS(J2+2, I2+3)=HESS(I2+3, J2+2)
194:          HESS(J2+3, I2+3)=HESS(I2+3, J2+3)175:          HESS(J2+3, I2+3)=HESS(I2+3, J2+3)
195: 176: 
196:       END DO177:       END DO
197:    END DO178:    END DO
198: 179: 
199: END SUBROUTINE180: END SUBROUTINE
200:  
201: SUBROUTINE STEALTHY_EIGEN_TEST(EN, GRAD) 
202:    USE COMMONS 
203:    USE KEY 
204:    USE MODHESS 
205:    IMPLICIT NONE 
206:    INTEGER I, NZEV, INFO, LSTECUT, HSTECUT 
207:    DOUBLE PRECISION GRAD(3*NATOMS), EIGEN(3*NATOMS), STTEMP(9*NATOMS) 
208:    DOUBLE PRECISION STRMS, EN 
209:  
210:    CALL DSYEV('N', 'U', 3*NATOMS, HESS, 3*NATOMS, EIGEN, STTEMP, 9*NATOMS, INFO) 
211:    STRMS=0 
212:  
213:    !write(*,*) EIGEN 
214:    DO I=1, NATOMS 
215:       STRMS=STRMS+GRAD(I)**2 
216:    ENDDO 
217:  
218:    STRMS=DSQRT(STRMS/NATOMS) 
219:    !write(*,'(T2,A,G20.10)') 'stealthy> Calculated RMS',STRMS 
220:  
221:    IF (STRMS.LE.GMAX) THEN 
222:       IF (EN.LT.1.0D-8) THEN 
223:          NZEV=NATOMS*3-STM*2 
224:          PRINT '(A)',' stealthy> Ground state' 
225:          !write(*,'(T2,A,I,G20.10)') 'stealthy> Test highest zero eigenvalue', NZEV, EIGEN(NZEV) 
226:          !write(*,'(T2,A,I,G20.10)') 'stealthy> Test lowest non-zero eigenvalue', NZEV+1, EIGEN(NZEV+1) 
227:          IF (EIGEN(1).LT.-1.0D-8) THEN 
228:             PRINT '(A,G20.10)',' stealthy> Error with lowest eigenvalue', EIGEN(1) 
229:          ENDIF 
230:          IF (EIGEN(NZEV+1).LT.1.0D-8) THEN 
231:             PRINT '(A,G20.10)',' stealthy> Error with lowest non-zero eigenvalue',  EIGEN(NZEV+1) 
232:          ENDIF 
233:          IF (EIGEN(NZEV).GT.1.0D-8) THEN 
234:             PRINT '(A,G20.10)',' stealthy> Error with highest zero eigenvalue', EIGEN(NZEV) 
235:          ENDIF 
236:       ELSE 
237:          NZEV=3 
238:          PRINT '(A)',' stealthy> Local minimum' 
239:          CALL ST_EIGEN_SORT(EIGEN, NZEV, LSTECUT, HSTECUT) 
240:          !write(*,'(T2,A,I,G20.10)') 'stealthy> Test highest zero eigenvalue', LSTECUT, EIGEN(LSTECUT) 
241:          !write(*,'(T2,A,I,G20.10)') 'stealthy> Test lowest non-zero eigenvalue', HSTECUT, EIGEN(HSTECUT) 
242:          IF (EIGEN(1).LT.-1.0D-8) THEN 
243:             PRINT '(A,G20.10)',' stealthy> Error with lowest eigenvalue', EIGEN(1) 
244:          ENDIF 
245:          IF (EIGEN(HSTECUT).LT.1.0D-8) THEN 
246:             PRINT '(A,G20.10)',' stealthy> Error with lowest non-zero eigenvalue', EIGEN(HSTECUT) 
247:          ENDIF 
248:          IF (EIGEN(LSTECUT).GT.1.0D-8) THEN 
249:             PRINT '(A,G20.10)',' stealthy> Error with highest zero eigenvalue', EIGEN(LSTECUT) 
250:          ENDIF 
251:       ENDIF 
252:    ELSE 
253:       NZEV=3 
254:       PRINT '(A)',' stealthy> Non-minimum configuration' 
255:       CALL ST_EIGEN_SORT(EIGEN, NZEV, LSTECUT, HSTECUT) 
256:       !write(*,'(T2,A,I,G20.10)') 'stealthy> Test highest zero eigenvalue', LSTECUT, EIGEN(LSTECUT) 
257:       !write(*,'(T2,A,I,G20.10)') 'stealthy> Test lowest non-zero eigenvalue', HSTECUT, EIGEN(HSTECUT+1) 
258:       IF (EIGEN(HSTECUT).LT.1.0D-8) THEN 
259:          PRINT '(A,G20.10)',' stealthy> Error with lowest non-zero eigenvalue', EIGEN(HSTECUT) 
260:       ENDIF 
261:       IF (EIGEN(LSTECUT).GT.1.0D-8) THEN 
262:          PRINT '(A,G20.10)',' stealthy> Error with highest zero eigenvalue', EIGEN(LSTECUT) 
263:       ENDIF 
264:    ENDIF 
265:  
266: END SUBROUTINE 
267:  
268: SUBROUTINE ST_EIGEN_SORT(EVARY, NLING, ELCUT, EHCUT) 
269:    USE COMMONS 
270:    IMPLICIT NONE 
271:    LOGICAL EQYN 
272:    INTEGER I, J, K, L, ELCUT, EHCUT, NLING 
273:    INTEGER,ALLOCATABLE :: BNNUM(:) 
274:    DOUBLE PRECISION EVARY(3*NATOMS), EVABS(3*NATOMS) 
275:    DOUBLE PRECISION MYMIN 
276:  
277:    ALLOCATE(BNNUM(NLING+1)) 
278:    EVABS=DABS(EVARY) 
279:    BNNUM=3*NATOMS 
280:    EQYN=.FALSE. 
281:     
282:    DO I=1, NLING+1 
283:       MYMIN=EVABS(3*NATOMS) 
284:       DO J=1, 3*NATOMS-1 
285:          DO K=1, I 
286:             IF (J==BNNUM(K)) EQYN=.TRUE. 
287:          ENDDO 
288:          IF (EQYN==.TRUE.) THEN 
289:             EQYN=.FALSE. 
290:             CYCLE 
291:          ENDIF 
292:          IF (EVABS(J).LT.MYMIN) THEN 
293:             MYMIN=EVABS(J) 
294:             BNNUM(I)=J 
295:          ENDIF 
296:       ENDDO 
297:       IF (I==NLING) ELCUT=BNNUM(I) 
298:       IF (I==NLING+1) EHCUT=BNNUM(I) 
299:    ENDDO 
300:  
301: END SUBROUTINE 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0