hdiff output

r30603/commons.f90 2016-06-16 09:30:15.357129417 +0100 r30602/commons.f90 2016-06-16 09:30:16.649146758 +0100
 32:      &        MYEUNIT, MYMUNIT, MYBUNIT, MYRUNIT, MYPUNIT, NFREEZETYPEA, & 32:      &        MYEUNIT, MYMUNIT, MYBUNIT, MYRUNIT, MYPUNIT, NFREEZETYPEA, &
 33:      &        TBPSTEPS, TBPCI, TBPBASIN, NTSITES, NRBGROUP, NZERO, PTMCDS_FRQ, PTMCDUMPENERFRQ, MONITORINT, NBLOCKS, & 33:      &        TBPSTEPS, TBPCI, TBPBASIN, NTSITES, NRBGROUP, NZERO, PTMCDS_FRQ, PTMCDUMPENERFRQ, MONITORINT, NBLOCKS, &
 34:      &        BINARY_EXAB_FRQ, NRESMIN, USERES, EXEQ, NONEDAPBC, STRUC, CHEMSHIFTITER, GRIDSIZE, MFETRUNS, BESTINVERT, GCNATOMS, & 34:      &        BINARY_EXAB_FRQ, NRESMIN, USERES, EXEQ, NONEDAPBC, STRUC, CHEMSHIFTITER, GRIDSIZE, MFETRUNS, BESTINVERT, GCNATOMS, &
 35:      &        GCINT, GCRELAX, MTARGETS, & 35:      &        GCINT, GCRELAX, MTARGETS, &
 36:      &        INTCONSEP, INTREPSEP, NCONSTRAINTON, CPREPSEP, CPCONSEP, NCONGEOM, & 36:      &        INTCONSEP, INTREPSEP, NCONSTRAINTON, CPREPSEP, CPCONSEP, NCONGEOM, &
 37:      &        NCPREPULSIVE, NCPCONSTRAINT, MAXCONUSE, INTCONSTEPS, INTRELSTEPS, INTSTEPS1, INTLJSTEPS, & 37:      &        NCPREPULSIVE, NCPCONSTRAINT, MAXCONUSE, INTCONSTEPS, INTRELSTEPS, INTSTEPS1, INTLJSTEPS, &
 38:      &        NTRAPPOW, MAXINTIMAGE, CHECKREPINTERVAL, INTFREEZEMIN, INTNTRIESMAX, INTIMAGEINCR, & 38:      &        NTRAPPOW, MAXINTIMAGE, CHECKREPINTERVAL, INTFREEZEMIN, INTNTRIESMAX, INTIMAGEINCR, &
 39:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, INTIMAGE, NREPULSIVE, & 39:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, INTIMAGE, NREPULSIVE, &
 40:      &        NNREPULSIVE, NCONSTRAINT, INTMUPDATE, DUMPINTEOSFREQ, DUMPINTXYZFREQ, & 40:      &        NNREPULSIVE, NCONSTRAINT, INTMUPDATE, DUMPINTEOSFREQ, DUMPINTXYZFREQ, &
 41:      &        LOCALPERMNEIGH, LOCALPERMMAXSEP, MAXNACTIVE, QCIPERMCHECKINT, & 41:      &        LOCALPERMNEIGH, LOCALPERMMAXSEP, MAXNACTIVE, QCIPERMCHECKINT, &
 42:      &        MLPIN, MLPSTART, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, DJWRBID, NHEXAMERS, QCIADDREP, QCIBONDS, QCISECOND 42:      &        MLPIN, MLPSTART, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, DJWRBID, NHEXAMERS
 43:   43:  
 44:       DOUBLE PRECISION RHO, GAMMA, SIG, SCEPS, SCC, TOLB, T12FAC, XMOVERENORM, RESIZE, QTSALLIS, & 44:       DOUBLE PRECISION RHO, GAMMA, SIG, SCEPS, SCC, TOLB, T12FAC, XMOVERENORM, RESIZE, QTSALLIS, &
 45:      &                 CQMAX, RADIUS, BQMAX,  MAXBFGS, DECAYPARAM, SYMTOL1, SYMTOL2, SYMTOL3, SYMTOL4, SYMTOL5, PGSYMTOLS(3),& 45:      &                 CQMAX, RADIUS, BQMAX,  MAXBFGS, DECAYPARAM, SYMTOL1, SYMTOL2, SYMTOL3, SYMTOL4, SYMTOL5, PGSYMTOLS(3),&
 46:      &                 ECONV, TOLD, TOLE, SYMREM(120,3,3), GMAX, CUTOFF, PCUT, EXPFAC, EXPD, CENTX, CENTY, CENTZ, & 46:      &                 ECONV, TOLD, TOLE, SYMREM(120,3,3), GMAX, CUTOFF, PCUT, EXPFAC, EXPD, CENTX, CENTY, CENTZ, &
 47:      &                 BOXLX, BOXLY, BOXLZ, BOX3D(3), PCUTOFF, SUPSTEP, SQUEEZER, SQUEEZED, COOPCUT, STOCKMU, STOCKLAMBDA, & 47:      &                 BOXLX, BOXLY, BOXLZ, BOX3D(3), PCUTOFF, SUPSTEP, SQUEEZER, SQUEEZED, COOPCUT, STOCKMU, STOCKLAMBDA, &
 48:      &                 TFAC(3), RMS, TEMPS, SACCRAT, CEIG, PNEWJUMP, EAMP, DISTFAC, ODDCHARGE, COULQ, COULSWAP, & 48:      &                 TFAC(3), RMS, TEMPS, SACCRAT, CEIG, PNEWJUMP, EAMP, DISTFAC, ODDCHARGE, COULQ, COULSWAP, &
 49:      &                 COULTEMP, APP, AMM, APM, XQP, XQM, ALPHAP, ALPHAM, ZSTAR, K_COMP, DGUESS, GUIDECUT, EFAC,&  49:      &                 COULTEMP, APP, AMM, APM, XQP, XQM, ALPHAP, ALPHAM, ZSTAR, K_COMP, DGUESS, GUIDECUT, EFAC,& 
 50:      &                 TRENORM, HISTMIN, HISTMAX, HISTFAC, EPSSPHERE, FINALCUTOFF, SHELLPROB, RINGROTSCALE, & 50:      &                 TRENORM, HISTMIN, HISTMAX, HISTFAC, EPSSPHERE, FINALCUTOFF, SHELLPROB, RINGROTSCALE, &
 51:      &                 HISTFACMUL, HPERCENT, AVOIDDIST, MAXERISE, MAXEFALL, TSTART, MATDIFF, STICKYSIG, SDTOL, & 51:      &                 HISTFACMUL, HPERCENT, AVOIDDIST, MAXERISE, MAXEFALL, TSTART, MATDIFF, STICKYSIG, SDTOL, &
 52:      &                 MinimalTemperature, MaximalTemperature, SwapProb, hdistconstraint, COREFRAC, TSTAR, & 52:      &                 MinimalTemperature, MaximalTemperature, SwapProb, hdistconstraint, COREFRAC, TSTAR, &
 75:      &                 PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2),PYA11(3),PYA21(3),PYA12(3),PYA22(3), & 75:      &                 PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2),PYA11(3),PYA21(3),PYA12(3),PYA22(3), &
 76:      &                 PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, PYCFTHRESH, LJSITECOORDS(3), & 76:      &                 PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, PYCFTHRESH, LJSITECOORDS(3), &
 77:      &                 MSTART,MFINISH,MBSTART1,MBFINISH1,MBSTART2,MBFINISH2,MBHEIGHT1,MBHEIGHT2,ME1,ME2,ME3, & 77:      &                 MSTART,MFINISH,MBSTART1,MBFINISH1,MBSTART2,MBFINISH2,MBHEIGHT1,MBHEIGHT2,ME1,ME2,ME3, &
 78:      &                 BSPTQMAX, BSPTQMIN, PFORCE, CSMNORM, CSMGUIDENORM, CSMEPS, PERCCUT, & 78:      &                 BSPTQMAX, BSPTQMIN, PFORCE, CSMNORM, CSMGUIDENORM, CSMEPS, PERCCUT, &
 79:      &                 LOWESTE, PERTSTEP, GCPLUS, & 79:      &                 LOWESTE, PERTSTEP, GCPLUS, &
 80:      &                 KINT, INTFREEZETOL, IMSEPMIN, IMSEPMAX, CONCUTABS, CONCUTFRAC, & 80:      &                 KINT, INTFREEZETOL, IMSEPMIN, IMSEPMAX, CONCUTABS, CONCUTFRAC, &
 81:      &                 LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, & 81:      &                 LPDGEOMDIFFTOL, INTCONFRAC, MAXCONE, INTRMSTOL, BFGSTSTOL, ORBITTOL, &
 82:      &                 INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, & 82:      &                 INTCONSTRAINTTOL, INTCONSTRAINTDEL, RBCUTOFF, INTCONSTRAINTREP, INTCONSTRAINREPCUT, &
 83:      &                 INTLJTOL, INTLJDEL, INTLJEPS, REPCON, INTDGUESS, CHECKREPCUTOFF, INTMINFAC, FREEZETOL, & 83:      &                 INTLJTOL, INTLJDEL, INTLJEPS, REPCON, INTDGUESS, CHECKREPCUTOFF, INTMINFAC, FREEZETOL, &
 84:      &                 LOCALPERMCUT, LOCALPERMCUT2, INTCONCUT, QCIRADSHIFT, MLPLAMBDA, & 84:      &                 LOCALPERMCUT, LOCALPERMCUT2, INTCONCUT, QCIRADSHIFT, MLPLAMBDA, &
 85:      &                 CAPSIDRHO,CAPSIDEPS,SIGMAPENT,RADPENT,SIGMAHEX,RADHEX,SIGMAPH, KLIM, SCA, & 85:      &                 CAPSIDRHO,CAPSIDEPS,SIGMAPENT,RADPENT,SIGMAHEX,RADHEX,SIGMAPH, KLIM, SCA
 86:      &                 QCIADDREPCUT, QCIADDREPEPS 
 87:  86: 
 88:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, & 87:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, &
 89:      &        OTPT, LJMFT, STRANDT, PAHT, SWT, MSTRANST, STOCKT, STICKYT, BLNT, MYSDT, FREEZERES, CENTXY, & 88:      &        OTPT, LJMFT, STRANDT, PAHT, SWT, MSTRANST, STOCKT, STICKYT, BLNT, MYSDT, FREEZERES, CENTXY, &
 90:      &        MSORIGT, SQUEEZET, PERIODIC, SCT, MSCT, MGUPTAT, RESIZET, TIP, RIGID, CALCQT, MPIT, JMT, LJCOULT, SETCENT, & 89:      &        MSORIGT, SQUEEZET, PERIODIC, SCT, MSCT, MGUPTAT, RESIZET, TIP, RIGID, CALCQT, MPIT, JMT, LJCOULT, SETCENT, &
 91:      &        SORTT, HIT, SAVEQ, PARALLELT, FIXD, RKMIN, BSMIN, PERMDIST, PERMOPT, BSWL, BSPT, BSPTRESTART, & 90:      &        SORTT, HIT, SAVEQ, PARALLELT, FIXD, RKMIN, BSMIN, PERMDIST, PERMOPT, BSWL, BSPT, BSPTRESTART, &
 92:      &        SYMMETRIZE, SYMMETRIZECSM, PRINT_PTGRP, PRINT_MINDATA, DUMPT, NEON, ARGON, P46, NORESET, TABOOT, EVSTEPT, PACHECO, DL_POLY, QUCENTRE, & 91:      &        SYMMETRIZE, SYMMETRIZECSM, PRINT_PTGRP, PRINT_MINDATA, DUMPT, NEON, ARGON, P46, NORESET, TABOOT, EVSTEPT, PACHECO, DL_POLY, QUCENTRE, &
 93:      &        STAR, PLUS, TWOPLUS, GROUND, DIPOLE, DFTBT, DFTBCT, SW, SUPERSTEP, EAMLJT, PBGLUET, TRACKDATAT, & 92:      &        STAR, PLUS, TWOPLUS, GROUND, DIPOLE, DFTBT, DFTBCT, SW, SUPERSTEP, EAMLJT, PBGLUET, TRACKDATAT, &
 94:      &        EAMALT, ALGLUET, MGGLUET, GUPTAT, LJATT, FST, DECAY, COOP, FIXBIN, GAUSST, QUENCHDOS, FIXDIHEFLAG, & 93:      &        EAMALT, ALGLUET, MGGLUET, GUPTAT, LJATT, FST, DECAY, COOP, FIXBIN, GAUSST, QUENCHDOS, FIXDIHEFLAG, &
 95:      &        FRAUSIT, ANGST, SELFT, STEPOUT, WENZEL, THRESHOLDT, THOMSONT, MULLERBROWNT, CHARMMENERGIES, & 94:      &        FRAUSIT, ANGST, SELFT, STEPOUT, WENZEL, THRESHOLDT, THOMSONT, MULLERBROWNT, CHARMMENERGIES, &
 96:      &        PROJ, RGCL2, TOSI, WELCH, AXTELL, AMBER, FIXIMAGE, BINARY, SHIFTCUT, ARNO, TUNNELT, TWOD, &  95:      &        PROJ, RGCL2, TOSI, WELCH, AXTELL, AMBER, FIXIMAGE, BINARY, SHIFTCUT, ARNO, TUNNELT, TWOD, & 
110:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&109:      &        HARMONICDONTMOVE, DUMPUNIQUE, FREEZESAVE, TBP, RBSYMT, PTMCDUMPSTRUCT, PTMCDUMPENERT, PYCOLDFUSION, MONITORT,&
111:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &110:      &        CHARMMDFTBT, PERMINVOPT, BLOCKMOVET, MAXERISE_SET, PYT, BINARY_EXAB, CHIROT, SANDBOXT, &
112:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &111:      &        RESERVOIRT, DISTOPT, ONEDAPBCT, ONEDPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, THREEDPBCT, RATIOT, &
113:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &112:      &        PTRANDOM, PTINTERVAL, PTSINGLE, PTSETS, CHEMSHIFT, CHEMSHIFT2, CSH, DEBUGss2029, UNIFORMMOVE, RANSEEDT, &
114:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &113:      &        TTM3T, NOINVERSION, RIGIDCONTOURT, UPDATERIGIDREFT, HYBRIDMINT, COMPRESSRIGIDT, MWFILMT, &
115:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &114:      &        SUPPRESST, MFETT, POLIRT, QUIPT, SWPOTT, MWPOTT, REPMATCHT, GLJT, MLJT, READMASST, SPECMASST, NEWTSALLIST, &
116:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 115:      &        PHI4MODELT, CUDAT, CUDATIMET, AMBER12T, ENERGY_DECOMPT, NEWMOVEST, DUMPMINT, MBPOLT, MOLECULART, GCBHT, SEMIGRAND_MUT, USEROT, & 
117:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &116:      &        SAVEMULTIMINONLY, GRADPROBLEMT, INTLJT, CONDATT, QCIPERMCHECK, &
118:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &117:      &        INTCONSTRAINTT, INTFREEZET, CHECKCONINT, CONCUTABST, CONCUTFRACT, INTERPCOSTFUNCTION, &
119:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &118:      &        RBAAT, FREEZENODEST, DUMPINTEOS, DUMPINTXYZ, QCIPOTT, QCIPOT2T, INTSPRINGACTIVET, LPERMDIST, LOCALPERMDIST, QCIRADSHIFTT, &
120:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT, QCINOREPINT119:      &        MLP3T, MKTRAPT, MLPB3T, MLPB3NEWT, MULTIPOTT, QCIAMBERT, MLPNEWREG, DJWRBT, STEALTHYT, LJADDT
121: !120: !
122:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 121:       DOUBLE PRECISION, ALLOCATABLE :: SEMIGRAND_MU(:) 
123:       DOUBLE PRECISION, ALLOCATABLE:: ATMASS(:)122:       DOUBLE PRECISION, ALLOCATABLE:: ATMASS(:)
124:       DOUBLE PRECISION, ALLOCATABLE:: SPECMASS(:) 123:       DOUBLE PRECISION, ALLOCATABLE:: SPECMASS(:) 
125: 124: 
126: ! csw34> FREEZEGROUP variables125: ! csw34> FREEZEGROUP variables
127: !126: !
128:       INTEGER :: GROUPCENTRE127:       INTEGER :: GROUPCENTRE
129:       DOUBLE PRECISION :: GROUPRADIUS128:       DOUBLE PRECISION :: GROUPRADIUS
130:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE129:       CHARACTER (LEN=2) :: FREEZEGROUPTYPE


r30603/congrad.f90 2016-06-16 09:30:15.561132153 +0100 r30602/congrad.f90 2016-06-16 09:30:16.869149713 +0100
   1: !   NEB module is an implementati on of the nudged elastic band method for performing double-ended pathway searches.
   2: !   Copyright (C) 2003-2006 Semen A. Trygubenko and David J. Wales
   3: !   This file is part of NEB module. NEB module is part of OPTIM.
   4: !
  1: !   OPTIM is free software; you can redistribute it and/or modify  5: !   OPTIM is free software; you can redistribute it and/or modify
  2: !   it under the terms of the GNU General Public License as published by  6: !   it under the terms of the GNU General Public License as published by
  3: !   the Free Software Foundation; either version 2 of the License, or  7: !   the Free Software Foundation; either version 2 of the License, or
  4: !   (at your option) any later version.  8: !   (at your option) any later version.
  5: !  9: !
  6: !   OPTIM is distributed in the hope that it will be useful, 10: !   OPTIM is distributed in the hope that it will be useful,
  7: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of
  8: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  9: !   GNU General Public License for more details. 13: !   GNU General Public License for more details.
 10: ! 14: !
502: 506: 
503: !507: !
504: ! This version of congrad tests for an internal minimum in the508: ! This version of congrad tests for an internal minimum in the
505: ! constraint distances as well as the repulsions.509: ! constraint distances as well as the repulsions.
506: !510: !
507: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)511: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
508: USE COMMONS, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &512: USE COMMONS, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
509:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &513:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
510:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &514:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &
511:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &515:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &
512:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, NATOMS, DEBUG, MYUNIT, INTMINFAC, INTSPRINGACTIVET516:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, NATOMS, DEBUG, MYUNIT, INTMINFAC, INTSPRINGACTIVET 
513: IMPLICIT NONE517: IMPLICIT NONE
514:            518:            
515: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2)519: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2)
516: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS520: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS
517: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1521: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
518: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)522: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
519: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI523: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
520: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)524: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)
521: LOGICAL NOINT, LPRINT525: LOGICAL NOINT, LPRINT
522: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2)526: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2)
897:    RMS=RMS+RMSIMAGE(J1)901:    RMS=RMS+RMSIMAGE(J1)
898:    IF (LPRINT) WRITE(MYUNIT, '(A,I6,2G20.10)') ' congrad2> J1,EEE,RMSIMAGE=',J1,EEE(J1),RMSIMAGE(J1)902:    IF (LPRINT) WRITE(MYUNIT, '(A,I6,2G20.10)') ' congrad2> J1,EEE,RMSIMAGE=',J1,EEE(J1),RMSIMAGE(J1)
899: ENDDO903: ENDDO
900: IF (INTIMAGE.NE.0) THEN904: IF (INTIMAGE.NE.0) THEN
901:    RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))905:    RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))
902: ENDIF906: ENDIF
903: !907: !
904: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,908: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points,
905: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)909: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2)
906: !910: !
907:  
908: ETOTAL=0.0D0911: ETOTAL=0.0D0
909: MAXINT=-1.0D100912: MAXINT=-1.0D100
910: MININT=1.0D100913: MININT=1.0D100
911: DO J1=2,INTIMAGE+1914: DO J1=2,INTIMAGE+1
912:    ETOTAL=ETOTAL+EEE(J1)915:    ETOTAL=ETOTAL+EEE(J1)
913: !  IF (DEBUG) PRINT '(A,I6,A,4G15.5)',' congrad2> con/rep and con/rep int image ', &916: !  IF (DEBUG) PRINT '(A,I6,A,4G15.5)',' congrad2> con/rep and con/rep int image ', &
914: ! &      J1,' ',CONE(J1),REPE(J1),CONEINT(J1),REPEINT(J1)917: ! &      J1,' ',CONE(J1),REPE(J1),CONEINT(J1),REPEINT(J1)
915:    IF (CONEINT(J1)+REPEINT(J1).LT.MININT) THEN918:    IF (CONEINT(J1)+REPEINT(J1).LT.MININT) THEN
916:       MININT=CONEINT(J1)+REPEINT(J1)919:       MININT=CONEINT(J1)+REPEINT(J1)
917:       NMININT=J1920:       NMININT=J1
966:    r1amr1bdr2amr2bsq=r1amr1bdr2amr2b**2969:    r1amr1bdr2amr2bsq=r1amr1bdr2amr2b**2
967:    r1amr1bsq=(r1ax - r1bx)**2 + (r1ay - r1by)**2 + (r1az - r1bz)**2970:    r1amr1bsq=(r1ax - r1bx)**2 + (r1ay - r1by)**2 + (r1az - r1bz)**2
968:    r2amr2bsq=(r2ax - r2bx)**2 + (r2ay - r2by)**2 + (r2az - r2bz)**2971:    r2amr2bsq=(r2ax - r2bx)**2 + (r2ay - r2by)**2 + (r2az - r2bz)**2
969:    DSQI=MAX((-r1amr1bdr2amr2bsq + r1amr1bsq*r2amr2bsq)/r1apr2bmr2amr1bsq,0.0D0)972:    DSQI=MAX((-r1amr1bdr2amr2bsq + r1amr1bsq*r2amr2bsq)/r1apr2bmr2amr1bsq,0.0D0)
970:    DINT=SQRT(DSQI)973:    DINT=SQRT(DSQI)
971: ENDIF974: ENDIF
972: 975: 
973: RETURN976: RETURN
974: 977: 
975: END SUBROUTINE INTMINONLY978: END SUBROUTINE INTMINONLY
976:  
977: ! 
978: ! Call this version for additional repulsive terms between constraints. 
979: ! 
980: SUBROUTINE CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
981: USE COMMONS, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, & 
982:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, & 
983:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,INTIMAGE, KINT, IMSEPMAX, ATOMACTIVE, QCINOREPINT, & 
984:   &            INTFREEZET, INTFROZEN, CONCUTLOCAL, CONCUT, CONCUTABST, CONCUTABS, CONCUTFRACT, CONCUTFRAC, & 
985:   &  NATOMS, DEBUG, MYUNIT, INTMINFAC, FREEZENODEST, INTSPRINGACTIVET, QCIADDREP, QCIADDREPCUT, QCIADDREPEPS, QCIBONDS 
986: USE PORFUNCS 
987: IMPLICIT NONE 
988:             
989: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,J3,J4,J5,NINTMIN,NINTMIN2 
990: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS 
991: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1 
992: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3) 
993: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI 
994: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2) 
995: LOGICAL NOINT 
996: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL 
997: LOGICAL IMGFREEZE(INTIMAGE) 
998: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3), X1, X2, Y1, Y2, Z1, Z2 
999:  
1000: EEE(1:INTIMAGE+2)=0.0D0 
1001: CONE(1:INTIMAGE+2)=0.0D0 
1002: REPE(1:INTIMAGE+2)=0.0D0 
1003: REPEINT(1:INTIMAGE+2)=0.0D0 
1004: NREPINT(1:INTIMAGE+2)=0 
1005: GGG(1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 
1006: ECON=0.0D0; EREP=0.0D0 
1007:  
1008: IF (QCIADDREP.LT.1) THEN 
1009:    WRITE(MYUNIT,'(A,I6)') 'congrad3> ERROR congrad3 called with no QCIADDREP=',QCIADDREP 
1010:    STOP 
1011: ENDIF 
1012: ! 
1013: !  Constraint energy and forces. 
1014: ! 
1015: OPEN(UNIT=852,FILE='test.xyz',STATUS='UNKNOWN') 
1016: INTCONST=QCIADDREPCUT**3 
1017: DO J2=1,NCONSTRAINT 
1018:    IF (.NOT.CONACTIVE(J2)) CYCLE 
1019:    CCLOCAL=CONCUTLOCAL(J2) 
1020:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS 
1021:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2) 
1022:    DO J1=2,INTIMAGE+1 
1023:       NI1=(3*NATOMS)*(J1-1)+3*(CONI(J2)-1) 
1024:       NJ1=(3*NATOMS)*(J1-1)+3*(CONJ(J2)-1) 
1025:       R2AX=XYZ(NI1+1); R2AY=XYZ(NI1+2); R2AZ=XYZ(NI1+3) 
1026:       R2BX=XYZ(NJ1+1); R2BY=XYZ(NJ1+2); R2BZ=XYZ(NJ1+3) 
1027:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2) 
1028:       IF (ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL) THEN  
1029:          DUMMY=D2-CONDISTREFLOCAL(J2)   
1030:          G2(1)=(R2AX-R2BX)/D2 
1031:          G2(2)=(R2AY-R2BY)/D2 
1032:          G2(3)=(R2AZ-R2BZ)/D2 
1033:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3) 
1034:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2) 
1035:          EEE(J1)=EEE(J1)  +DUMMY 
1036:          ECON=ECON        +DUMMY 
1037:          CONE(J1)=CONE(J1)+DUMMY 
1038:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3) 
1039:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3) 
1040:       ENDIF 
1041:       IF (J2.GT.QCIBONDS) CYCLE 
1042:       DO J3=J2+1,QCIBONDS 
1043:          IF (.NOT.CONACTIVE(J3)) CYCLE 
1044:          IF (CONI(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom 
1045:          IF (CONI(J3).EQ.CONJ(J2)) CYCLE ! no extra terms for bonds with a common atom 
1046:          IF (CONJ(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom 
1047:          IF (CONJ(J3).EQ.CONJ(J2)) CYCLE ! no extra terms for bonds with a common atom 
1048:          NI2=(3*NATOMS)*(J1-1)+3*(CONI(J3)-1) 
1049:          NJ2=(3*NATOMS)*(J1-1)+3*(CONJ(J3)-1) 
1050:          DO J4=1,QCIADDREP 
1051:             X1=(J4*XYZ(NI1+1)+(QCIADDREP+1-J4)*XYZ(NJ1+1))/(QCIADDREP+1.0D0) 
1052:             Y1=(J4*XYZ(NI1+2)+(QCIADDREP+1-J4)*XYZ(NJ1+2))/(QCIADDREP+1.0D0) 
1053:             Z1=(J4*XYZ(NI1+3)+(QCIADDREP+1-J4)*XYZ(NJ1+3))/(QCIADDREP+1.0D0) 
1054:             DO J5=1,QCIADDREP 
1055:                X2=(J5*XYZ(NI2+1)+(QCIADDREP+1-J5)*XYZ(NJ2+1))/(QCIADDREP+1.0D0) 
1056:                Y2=(J5*XYZ(NI2+2)+(QCIADDREP+1-J5)*XYZ(NJ2+2))/(QCIADDREP+1.0D0) 
1057:                Z2=(J5*XYZ(NI2+3)+(QCIADDREP+1-J5)*XYZ(NJ2+3))/(QCIADDREP+1.0D0) 
1058:                D2=SQRT((X1-X2)**2+(Y1-Y2)**2+(Z1-Z2)**2) 
1059:                IF (D2.LT.QCIADDREPCUT) THEN  
1060:                   D12=D2**2 
1061:                   DUMMY=QCIADDREPEPS*(1.0D0/D12+(2.0D0*D2-3.0D0*QCIADDREPCUT)/INTCONST) 
1062:                   EEE(J1)=EEE(J1)+DUMMY 
1063:                   REPE(J1)=REPE(J1)+DUMMY 
1064:                   EREP=EREP+DUMMY 
1065: !                 WRITE(MYUNIT,'(A,4I6,A,2I6,A,2G20.10)') 'congrad3> atoms: ',CONI(J2),CONJ(J2),CONI(J3),CONJ(J3), & 
1066: ! &                     ' sites ',J4,J5,' dist,erep=',D2,DUMMY    
1067:                   DUMMY=-2.0D0*QCIADDREPEPS*(1.0D0/(D2*D12)-1.0D0/INTCONST) 
1068:                   G2(1)=(X1-X2)/D2 
1069:                   G2(2)=(Y1-Y2)/D2 
1070:                   G2(3)=(Z1-Z2)/D2 
1071:                   REPGRAD(1:3)=DUMMY*G2(1:3) 
1072:                   GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+J4*REPGRAD(1:3)/(QCIADDREP+1.0D0) ! forces on the four atoms involved in the two constraints 
1073:                   GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)+(QCIADDREP+1-J4)*REPGRAD(1:3)/(QCIADDREP+1.0D0) 
1074:                   GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)-J5*REPGRAD(1:3)/(QCIADDREP+1.0D0) 
1075:                   GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-(QCIADDREP+1-J5)*REPGRAD(1:3)/(QCIADDREP+1.0D0) 
1076:                ENDIF 
1077:             ENDDO 
1078:          ENDDO 
1079:       ENDDO 
1080:    ENDDO 
1081: ENDDO 
1082: CLOSE(852) 
1083:  
1084: GGG(1:(3*NATOMS))=0.0D0                                      ! can delete when loop range above changes 
1085: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes 
1086:  
1087: ! INTCONST=INTCONSTRAINREPCUT**13 
1088:  
1089: DO J2=1,NNREPULSIVE 
1090: !  INTCONST=NREPCUT(J2)**13 
1091:    INTCONST=NREPCUT(J2)**3 
1092:    DO J1=2,INTIMAGE+2 
1093: !  DO J1=1,INTIMAGE+2 ! can change when zero energies are confirmed for end images 
1094:       IF (FREEZENODEST) THEN 
1095:          IF (J1.EQ.2) THEN 
1096:             IF (IMGFREEZE(1)) CYCLE 
1097:          ELSE IF (J1.EQ.INTIMAGE+2) THEN 
1098:             IF (IMGFREEZE(INTIMAGE)) CYCLE 
1099:          ELSE 
1100:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE 
1101:          ENDIF 
1102:       ENDIF 
1103:       IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN 
1104:          WRITE(MYUNIT, '(A,I6,A,2I6)') ' congrad> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2) 
1105:          STOP 
1106:       ENDIF 
1107: !     WRITE(MYUNIT,'(A,2I8,6G20.10)') 'congrad2> B J1,J2,GGG(1:6)=',J1,J2,GGG(1:6) 
1108:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1) 
1109:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1) 
1110:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3) 
1111:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3) 
1112:       D2=SQRT((R2AX-R2BX)**2+(R2AY-R2BY)**2+(R2AZ-R2BZ)**2) 
1113:       IF (D2.LT.NREPCUT(J2)) THEN ! term for image J1 
1114: !        D12=D2**12 
1115:          D12=D2**2 
1116: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST) 
1117:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST) 
1118:          EEE(J1)=EEE(J1)+DUMMY 
1119:          REPE(J1)=REPE(J1)+DUMMY 
1120:          EREP=EREP+DUMMY 
1121: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST) 
1122:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST) 
1123:          G2(1)=(R2AX-R2BX)/D2 
1124:          G2(2)=(R2AY-R2BY)/D2 
1125:          G2(3)=(R2AZ-R2BZ)/D2 
1126:          REPGRAD(1:3)=DUMMY*G2(1:3) 
1127:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3) 
1128:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3) 
1129:       ENDIF 
1130: ! 
1131: ! For internal minima we are counting edges.  
1132: ! Edge J1 is between images J1-1 and J1, starting from J1=2. 
1133: ! Energy contributions are shared evenly, except for 
1134: ! edge 1, which is assigned to image 2, and edge INTIMAGE+1, which 
1135: ! is assigned to image INTIMAGE+1. Gradients are set to zero for 
1136: ! the end images. 
1137: ! 
1138:       IF (J1.EQ.1) CYCLE 
1139:       IF (QCINOREPINT) CYCLE 
1140:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1) 
1141:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1) 
1142:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3) 
1143:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3) 
1144: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN 
1145:       IF ((r1ax-r1bx-r2ax+r2bx)**2+(r1ay-r1by-r2ay+r2by)**2+(r1az-r1bz-r2az+r2bz)**2.LT.1.0D-10) THEN 
1146: !        WRITE(MYUNIT, '(A,I6,A,2I6)') 'B repulsion number ',J2, ' between ',NREPI(J2),NREPJ(J2) 
1147: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ=',R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ 
1148: !        WRITE(MYUNIT, '(A,6F15.10)') 'R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ=',R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ 
1149: !        WRITE(MYUNIT,'(A,7I10)') 'congrad> J2,NI1,NJ1,NI2,NJ2,NREPI,NREPJ=',J2,NI1,NJ1,NI2,NJ2,NREPI(J2),NREPJ(J2) 
1150: !        WRITE(MYUNIT,'(A,7I10)') 'frames ',J1-1,J1 
1151:          NOINT=.TRUE. 
1152:       ELSE 
1153:          CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, & 
1154:   &                 D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2)) 
1155:          IF (.NOT.NOINT) THEN 
1156: !           WRITE(MYUNIT,'(A,I6,A,I6,A,2I6,A,2G20.10)') 'congrad> internal minimum images ',J1-1,' and ',J1,' atoms: ',NREPI(J2),NREPJ(J2), & 
1157: ! &                        ' distance,cutoff=',DINT,NREPCUT(J2) 
1158:             NINTMIN=NINTMIN+1 
1159:          ENDIF 
1160:       ENDIF 
1161:       IF ((.NOT.NOINT).AND.(DINT.LT.NREPCUT(J2))) THEN 
1162:          NINTMIN2=NINTMIN2+1 
1163: !        D12=DSQI**6 
1164:          D12=DSQI 
1165: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST) 
1166:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST) 
1167:          IF (J1.EQ.2) THEN 
1168:             EEE(J1)=EEE(J1)+DUMMY 
1169:             REPEINT(J1)=REPEINT(J1)+DUMMY 
1170:             NREPINT(J1)=NREPINT(J1)+1 
1171:          ELSE IF (J1.LT.INTIMAGE+2) THEN 
1172:             EEE(J1)=EEE(J1)+DUMMY/2.0D0 
1173:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0 
1174:             REPEINT(J1)=REPEINT(J1)+DUMMY/2.0D0 
1175:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0 
1176:             NREPINT(J1)=NREPINT(J1)+1 
1177:             NREPINT(J1-1)=NREPINT(J1-1)+1 
1178:          ELSE IF (J1.EQ.INTIMAGE+2) THEN 
1179:             EEE(J1-1)=EEE(J1-1)+DUMMY 
1180:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY 
1181:             NREPINT(J1-1)=NREPINT(J1-1)+1 
1182:          ENDIF 
1183:          EREP=EREP+DUMMY 
1184: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST) 
1185:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(DINT*D12)-1.0D0/INTCONST) 
1186:          REPGRAD(1:3)=INTMINFAC*DUMMY*G1INT(1:3) 
1187: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), & 
1188: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2) 
1189: ! 
1190: ! Gradient contributions for image J1-1 
1191: ! 
1192:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3) 
1193:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3) 
1194:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3) 
1195: !        PRINT '(A,4I6,2G15.5)','in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), & 
1196: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2) 
1197: ! 
1198: ! Gradient contributions for image J1 
1199: ! 
1200:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3) 
1201:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3) 
1202:       ENDIF 
1203:    ENDDO 
1204: ENDDO 
1205: ! 
1206: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise 
1207: ! energy terms between images except for the end points. 
1208: ! 
1209: ESPRING=0.0D0 
1210: IF (KINT.NE.0.0D0) THEN 
1211:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1 
1212:       NI1=(3*NATOMS)*(J1-1) 
1213:       NI2=(3*NATOMS)*J1 
1214: ! 
1215: !  Edge between J1 and J1+1 
1216: ! 
1217:       DPLUS=0.0D0 
1218:       DO J2=1,NATOMS 
1219:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN  
1220:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 & 
1221:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 & 
1222:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2 
1223:          ENDIF 
1224:       ENDDO 
1225:       DPLUS=SQRT(DPLUS) 
1226:       IF (DPLUS.GT.IMSEPMAX) THEN 
1227: !        DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2 
1228:          DUMMY=KINT*0.5D0*DPLUS**2 
1229:          ESPRING=ESPRING+DUMMY 
1230: !        DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS 
1231:          DUMMY=KINT 
1232:          DO J2=1,NATOMS 
1233:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN  
1234:                SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)) 
1235:                GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3) 
1236:                GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3) 
1237:             ENDIF 
1238:          ENDDO 
1239:       ENDIF 
1240:    ENDDO 
1241: ENDIF 
1242: ! 
1243: ! Set gradients on frozen atoms to zero. 
1244: ! 
1245: IF (FREEZE) THEN 
1246:    DO J1=2,INTIMAGE+1   
1247:       DO J2=1,NATOMS 
1248:          IF (FROZEN(J2)) THEN 
1249:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+1)=0.0D0 
1250:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+2)=0.0D0 
1251:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+3)=0.0D0 
1252:          ENDIF 
1253:       ENDDO 
1254:    ENDDO 
1255: ENDIF 
1256: ! 
1257: ! Set gradients on locally frozen atoms to zero. 
1258: ! 
1259: IF (INTFREEZET) THEN 
1260:    DO J1=2,INTIMAGE+1   
1261:       DO J2=1,NATOMS 
1262:          IF (INTFROZEN(J2)) THEN 
1263:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+1)=0.0D0 
1264:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+2)=0.0D0 
1265:             GGG((3*NATOMS)*(J1-1)+3*(J2-1)+3)=0.0D0 
1266:          ENDIF 
1267:       ENDDO 
1268:    ENDDO 
1269: ENDIF 
1270: ! 
1271: ! Set gradients to zero for start and finish images. 
1272: ! 
1273: GGG(1:(3*NATOMS))=0.0D0 
1274: GGG((INTIMAGE+1)*(3*NATOMS)+1:(INTIMAGE+2)*(3*NATOMS))=0.0D0 
1275: RMS=0.0D0 
1276: DO J1=2,INTIMAGE+1 
1277:    RMSIM(J1)=0.0D0 
1278:    DO J2=1,(3*NATOMS) 
1279:       RMS=RMS+GGG((3*NATOMS)*(J1-1)+J2)**2 
1280:       RMSIM(J1)=RMSIM(J1)+GGG((3*NATOMS)*(J1-1)+J2)**2 
1281:    ENDDO 
1282:    RMSIM(J1)=SQRT(RMSIM(J1)/(3*NATOMS)) 
1283: ENDDO 
1284: IF (INTIMAGE.NE.0) THEN 
1285:    RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE)) 
1286: ENDIF 
1287: ! 
1288: ! For INTIMAGE images there are INTIMAGE+2 replicas including the end points, 
1289: ! and INTIMAGE+1 line segements, with associated energies stored in EEE(2:INTIMAGE+2) 
1290: ! 
1291: ETOTAL=0.0D0 
1292: MAXINT=-1.0D100 
1293: MININT=1.0D100 
1294: DO J1=2,INTIMAGE+1 
1295:    ETOTAL=ETOTAL+EEE(J1) 
1296: !  WRITE(MYUNIT, '(A,I6,A,3G20.10)') ' congrad> con/rep/RMS image ',J1,' ',CONE(J1),REPE(J1),RMSIM(J1) 
1297:    IF (REPEINT(J1).LT.MININT) THEN 
1298:       MININT=REPEINT(J1) 
1299:       NMININT=J1 
1300:    ENDIF 
1301:    IF (REPE(J1).GT.MAXINT) THEN 
1302:       MAXINT=REPE(J1) 
1303:       NMAXINT=J1 
1304:    ENDIF 
1305: ENDDO 
1306:  
1307: END SUBROUTINE CONGRAD3 


r30603/intlbfgs.f90 2016-06-16 09:30:15.781135129 +0100 r30602/intlbfgs.f90 2016-06-16 09:30:17.137153320 +0100
 21:      & INTRMSTOL, INTIMAGE, NREPMAX, NREPULSIVE, MUPDATE, INTDGUESS, & 21:      & INTRMSTOL, INTIMAGE, NREPMAX, NREPULSIVE, MUPDATE, INTDGUESS, &
 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, & 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, &
 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, & 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, &
 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, & 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, &
 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, & 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, &
 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, & 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, &
 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, & 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, &
 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, & 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, &
 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, & 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, &
 30:      & CONCUT, CONCUTLOCAL, NATOMS, DEBUG, STEP, MCSTEPS, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, & 30:      & CONCUT, CONCUTLOCAL, NATOMS, DEBUG, STEP, MCSTEPS, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, &
 31:      & NCONSTRAINTFIX, CONIFIX, CONJFIX, QCIPERMCHECK, QCIPERMCHECKINT, PERIODIC, TWOD, RIGID, BOXLX, BOXLY, BOXLZ, & 31:      & NCONSTRAINTFIX, CONIFIX, CONJFIX, QCIPERMCHECK, QCIPERMCHECKINT, PERIODIC, TWOD, RIGID, BOXLX, BOXLY, BOXLZ
 32:      & QCIADDREP 
 33:  32: 
 34: IMPLICIT NONE  33: IMPLICIT NONE 
 35:  34: 
 36: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points 35: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points
 37: INTEGER D, U 36: INTEGER D, U
 38: DOUBLE PRECISION DIST, DIST2, RMAT(3,3) 37: DOUBLE PRECISION DIST, DIST2, RMAT(3,3)
 39: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ 38: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ
 40: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE, JA1, JA2 39: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE, JA1, JA2
 41: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS), LDEBUG 40: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS), LDEBUG
 42: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 41: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 65:   &              EEETMP(:), MYGTMP(:), EEE(:), STEPIMAGE(:), & 64:   &              EEETMP(:), MYGTMP(:), EEE(:), STEPIMAGE(:), &
 66:   &              GTMP(:), DIAG(:), STP(:), SEARCHSTEP(:,:), GDIF(:,:), GLAST(:), XSAVE(:) 65:   &              GTMP(:), DIAG(:), STP(:), SEARCHSTEP(:,:), GDIF(:,:), GLAST(:), XSAVE(:)
 67: DOUBLE PRECISION, ALLOCATABLE, TARGET :: XYZ(:), GGG(:), DPTMP(:), D2TMP(:,:) 66: DOUBLE PRECISION, ALLOCATABLE, TARGET :: XYZ(:), GGG(:), DPTMP(:), D2TMP(:,:)
 68: ! saved interpolation 67: ! saved interpolation
 69: DOUBLE PRECISION, ALLOCATABLE :: BESTXYZ(:), BESTEEE(:) 68: DOUBLE PRECISION, ALLOCATABLE :: BESTXYZ(:), BESTEEE(:)
 70: INTEGER BESTINTIMAGE, NSTEPS, NITERUSE 69: INTEGER BESTINTIMAGE, NSTEPS, NITERUSE
 71: LOGICAL, ALLOCATABLE :: CHECKG(:), IMGFREEZE(:) 70: LOGICAL, ALLOCATABLE :: CHECKG(:), IMGFREEZE(:)
 72: LOGICAL READIMAGET 71: LOGICAL READIMAGET
 73: INTEGER LUNIT, GETUNIT 72: INTEGER LUNIT, GETUNIT
 74: CHARACTER(LEN=2) SDUMMY 73: CHARACTER(LEN=2) SDUMMY
 75: INTEGER JMAXEEE,JMAXRMS 
 76: DOUBLE PRECISION MAXEEE,MAXRMS 
 77:  74: 
 78: READIMAGET=.TRUE. 75: READIMAGET=.TRUE.
 79: READIMAGET=.FALSE. 76: READIMAGET=.FALSE.
 80: IF (READIMAGET) THEN 77: IF (READIMAGET) THEN
 81:    LUNIT=GETUNIT() 78:    LUNIT=GETUNIT()
 82:    OPEN(UNIT=LUNIT,FILE='restart.xyz',STATUS='OLD') 79:    OPEN(UNIT=LUNIT,FILE='restart.xyz',STATUS='OLD')
 83:    INTIMAGE=0 80:    INTIMAGE=0
 84: 653 CONTINUE 81: 653 CONTINUE
 85:    READ(LUNIT,*,END=654) NDUMMY 82:    READ(LUNIT,*,END=654) NDUMMY
 86:    READ(LUNIT,*)  83:    READ(LUNIT,*) 
451: NSTEPSMAX=INTSTEPS1448: NSTEPSMAX=INTSTEPS1
452: !449: !
453: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.450: ! Don;t want to redistribute images before even taking a step, so don;t call CHECKSEP.
454: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!451: ! Must call CHECKREP to initialise NNREULSIVE, NREPI, NREPJ, etc. SEGV otherwise on second cycle!
455: !452: !
456: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.453: ! To take BH-type steps in the QCI space, jump back here. Leave SWITCHED true.
457: !454: !
458: BESTWORST=1.0D100455: BESTWORST=1.0D100
459: 9876 CONTINUE456: 9876 CONTINUE
460: CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)457: CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
461: IF (QCIADDREP.GT.0) THEN458: IF (CHECKCONINT) THEN
462:    CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
463: ELSEIF (CHECKCONINT) THEN 
464:    CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)459:    CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
465: ELSE460: ELSE
466:    CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)461:    CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
467: ENDIF462: ENDIF
468: EOLD=ETOTAL463: EOLD=ETOTAL
469: GLAST(1:D)=G(1:D)464: GLAST(1:D)=G(1:D)
470: XSAVE(1:D)=X(1:D)465: XSAVE(1:D)=X(1:D)
471: 466: 
472: IF (ETOTAL/INTIMAGE.LT.COLDFUSIONLIMIT) THEN467: IF (ETOTAL/INTIMAGE.LT.COLDFUSIONLIMIT) THEN
473:    WRITE(MYUNIT,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=', &468:    WRITE(MYUNIT,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=', &
753:             CHECKG(1:(3*NATOMS)*(INTIMAGE+1))=.FALSE.748:             CHECKG(1:(3*NATOMS)*(INTIMAGE+1))=.FALSE.
754:             IMGFREEZE(1:INTIMAGE+1)=.FALSE.749:             IMGFREEZE(1:INTIMAGE+1)=.FALSE.
755:             EEE(1:INTIMAGE+3)=0.0D0750:             EEE(1:INTIMAGE+3)=0.0D0
756:             STEPIMAGE(1:INTIMAGE+1)=0.0D0751:             STEPIMAGE(1:INTIMAGE+1)=0.0D0
757: 752: 
758:             X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+2))753:             X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+2))
759:             G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+2))754:             G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+2))
760:             INTIMAGE=INTIMAGE+1755:             INTIMAGE=INTIMAGE+1
761:             D=(3*NATOMS)*INTIMAGE756:             D=(3*NATOMS)*INTIMAGE
762:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)757:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
763:             IF (QCIADDREP.GT.0) THEN758:             IF (CHECKCONINT) THEN
764:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
765:             ELSEIF (CHECKCONINT) THEN 
766:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)759:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
767:             ELSE760:             ELSE
768:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)761:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
769:             ENDIF762:             ENDIF
770: !           GOTO 864763: !           GOTO 864
771:          ENDIF764:          ENDIF
772:          IF ((DMIN.LT.IMSEPMIN).AND.(INTIMAGE.GT.1)) THEN765:          IF ((DMIN.LT.IMSEPMIN).AND.(INTIMAGE.GT.1)) THEN
773:             IF (JMIN.EQ.1) JMIN=2766:             IF (JMIN.EQ.1) JMIN=2
774:             WRITE(MYUNIT,'(A,I6,A,I6)') ' intlbfgs> Remove image ',JMIN767:             WRITE(MYUNIT,'(A,I6,A,I6)') ' intlbfgs> Remove image ',JMIN
775:             NITERUSE=0768:             NITERUSE=0
832:             CHECKG(1:(3*NATOMS)*(INTIMAGE-1))=.FALSE.825:             CHECKG(1:(3*NATOMS)*(INTIMAGE-1))=.FALSE.
833:             IMGFREEZE(1:INTIMAGE-1)=.FALSE.826:             IMGFREEZE(1:INTIMAGE-1)=.FALSE.
834:             EEE(1:INTIMAGE+1)=0.0D0827:             EEE(1:INTIMAGE+1)=0.0D0
835:             STEPIMAGE(1:INTIMAGE-1)=0.0D0828:             STEPIMAGE(1:INTIMAGE-1)=0.0D0
836: 829: 
837:             X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE))830:             X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE))
838:             G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE))831:             G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE))
839:             INTIMAGE=INTIMAGE-1832:             INTIMAGE=INTIMAGE-1
840:             D=(3*NATOMS)*INTIMAGE833:             D=(3*NATOMS)*INTIMAGE
841:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)834:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
842:             IF (QCIADDREP.GT.0) THEN835:             IF (CHECKCONINT) THEN
843:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
844:             ELSEIF (CHECKCONINT) THEN 
845:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)836:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
846:             ELSE837:             ELSE
847:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)838:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
848:             ENDIF839:             ENDIF
849: !           GOTO 864840: !           GOTO 864
850:          ENDIF841:          ENDIF
851:       ELSE842:       ELSE
852:          DMAX=-1.0D0843:          DMAX=-1.0D0
853:          ADMAX=-1.0D0844:          ADMAX=-1.0D0
854:          DMIN=HUGE(1.0D0)845:          DMIN=HUGE(1.0D0)
886: !            IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,I6,A,G20.10)')' intlbfgs> distance between images ', &877: !            IF (DEBUG) WRITE(MYUNIT,'(A,I6,A,I6,A,G20.10)')' intlbfgs> distance between images ', &
887: !  &                                                  J1,' and ',J1+1,' is ',DUMMY878: !  &                                                  J1,' and ',J1+1,' is ',DUMMY
888: !!           IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6)')' intlbfgs> largest atomic distance between images so far is ', &879: !!           IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6)')' intlbfgs> largest atomic distance between images so far is ', &
889: !! &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1880: !! &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1
890:          ENDDO881:          ENDDO
891:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &882:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &
892:   &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1,' total images=',INTIMAGE883:   &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1,' total images=',INTIMAGE
893:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &884:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &
894:   &                                                  DMAX,' for images ',JMAX,JMAX+1885:   &                                                  DMAX,' for images ',JMAX,JMAX+1
895:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,G20.10)') 'intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)886:          IF (DEBUG) WRITE(MYUNIT,'(A,G20.10,A,G20.10)') 'intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)
896: !        IF (SQRT(ADMAX).GT.IMSEPMAX) THEN887:          IF (SQRT(ADMAX).GT.IMSEPMAX) THEN
897: !           KINT=MIN(1.0D6,KINT*1.1D0)888:             KINT=MIN(1.0D6,KINT*1.1D0)
898: !        ELSE889:          ELSE
899: !           KINT=MAX(1.0D-6,KINT/1.1D0)890:             KINT=MAX(1.0D-6,KINT/1.1D0)
900: !        ENDIF891:          ENDIF
901:          WRITE(MYUNIT,'(A,G20.10)') 'intlbfgs> Spring constant is now ',KINT892:          WRITE(MYUNIT,'(A,G20.10)') 'intlbfgs> Spring constant is now ',KINT
902:       ENDIF893:       ENDIF
903:    ENDIF894:    ENDIF
904: !895: !
905: ! End of add/subtract images block.896: ! End of add/subtract images block.
906: !897: !
907:    IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN898:    IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN
908:       LDEBUG=.FALSE.899:       LDEBUG=.FALSE.
909:       DO J2=2,INTIMAGE+2900:       DO J2=2,INTIMAGE+2
910:          CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &901:          CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &
917: !     IF (INTIMAGE.GT.300) THEN908: !     IF (INTIMAGE.GT.300) THEN
918: !        WRITE(MYUNIT,'(A,2L5)') 'atom 375 intfrozen and atomactive: ',INTFROZEN(375),ATOMACTIVE(375)909: !        WRITE(MYUNIT,'(A,2L5)') 'atom 375 intfrozen and atomactive: ',INTFROZEN(375),ATOMACTIVE(375)
919: !        WRITE(MYUNIT,'(A,2L5)') 'atom 384 intfrozen and atomactive: ',INTFROZEN(384),ATOMACTIVE(384)910: !        WRITE(MYUNIT,'(A,2L5)') 'atom 384 intfrozen and atomactive: ',INTFROZEN(384),ATOMACTIVE(384)
920: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 284:',XYZ(3*400*283+3*374+1:3*400*283+3*374+3)911: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 284:',XYZ(3*400*283+3*374+1:3*400*283+3*374+3)
921: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 285:',XYZ(3*400*284+3*374+1:3*400*284+3*374+3)912: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 285:',XYZ(3*400*284+3*374+1:3*400*284+3*374+3)
922: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 286:',XYZ(3*400*285+3*374+1:3*400*285+3*374+3)913: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 286:',XYZ(3*400*285+3*374+1:3*400*285+3*374+3)
923: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 284:',XYZ(3*400*283+3*383+1:3*400*283+3*383+3)914: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 284:',XYZ(3*400*283+3*383+1:3*400*283+3*383+3)
924: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 285:',XYZ(3*400*284+3*383+1:3*400*284+3*383+3)915: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 285:',XYZ(3*400*284+3*383+1:3*400*284+3*383+3)
925: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 286:',XYZ(3*400*285+3*383+1:3*400*285+3*383+3)916: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 286:',XYZ(3*400*285+3*383+1:3*400*285+3*383+3)
926: !     ENDIF917: !     ENDIF
927:       IF (QCIADDREP.GT.0) THEN918:       IF (CHECKCONINT) THEN
928:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
929:       ELSEIF (CHECKCONINT) THEN 
930:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)919:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
931:       ELSE920:       ELSE
932:          CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)921:          CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
933:       ENDIF922:       ENDIF
934: !     IF (INTIMAGE.GT.300) THEN923: !     IF (INTIMAGE.GT.300) THEN
935: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 284 GGG:',GGG(3*400*283+3*374+1:3*400*283+3*374+3)924: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 284 GGG:',GGG(3*400*283+3*374+1:3*400*283+3*374+3)
936: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 285 GGG:',GGG(3*400*284+3*374+1:3*400*284+3*374+3)925: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 285 GGG:',GGG(3*400*284+3*374+1:3*400*284+3*374+3)
937: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 286 GGG:',GGG(3*400*285+3*374+1:3*400*285+3*374+3)926: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 375 in image 286 GGG:',GGG(3*400*285+3*374+1:3*400*285+3*374+3)
938: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 284 GGG:',GGG(3*400*283+3*383+1:3*400*283+3*383+3)927: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 284 GGG:',GGG(3*400*283+3*383+1:3*400*283+3*383+3)
939: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 285 GGG:',GGG(3*400*284+3*383+1:3*400*284+3*383+3)928: !        WRITE(MYUNIT,'(A,3G20.10)') 'atom 384 in image 285 GGG:',GGG(3*400*284+3*383+1:3*400*284+3*383+3)
968:          DO J4=2,INTIMAGE+1957:          DO J4=2,INTIMAGE+1
969:             IF (CHRMMT) CALL UPDATENBONDS(XYZ((3*NATOMS)*(J4-1)+1:(3*NATOMS)*J4))958:             IF (CHRMMT) CALL UPDATENBONDS(XYZ((3*NATOMS)*(J4-1)+1:(3*NATOMS)*J4))
970:             CALL POTENTIAL(XYZ((3*NATOMS)*(J4-1)+1:(3*NATOMS)*J4),GGG((3*NATOMS)*(J4-1)+1:(3*NATOMS)*J4),EEE(J4), &959:             CALL POTENTIAL(XYZ((3*NATOMS)*(J4-1)+1:(3*NATOMS)*J4),GGG((3*NATOMS)*(J4-1)+1:(3*NATOMS)*J4),EEE(J4), &
971:   &                                    .TRUE.,.FALSE.)960:   &                                    .TRUE.,.FALSE.)
972:             ETOTALTMP=ETOTALTMP+EEE(J4)961:             ETOTALTMP=ETOTALTMP+EEE(J4)
973:          ENDDO962:          ENDDO
974:       ENDIF963:       ENDIF
975:       EEETMP(1:INTIMAGE+2)=EEE(1:INTIMAGE+2)964:       EEETMP(1:INTIMAGE+2)=EEE(1:INTIMAGE+2)
976:       MYGTMP(1:D)=G(1:D)965:       MYGTMP(1:D)=G(1:D)
977:       IF (USEFRAC.LT.1.0D0) THEN966:       IF (USEFRAC.LT.1.0D0) THEN
978:          IF (QCIADDREP.GT.0) THEN967:          IF (CHECKCONINT) THEN
979:             CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
980:          ELSEIF (CHECKCONINT) THEN 
981:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)968:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
982:          ELSE969:          ELSE
983:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)970:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
984:          ENDIF971:          ENDIF
985:       ELSE972:       ELSE
986:          ETOTAL=0.0D0973:          ETOTAL=0.0D0
987:          G(1:D)=0.0D0974:          G(1:D)=0.0D0
988:       ENDIF975:       ENDIF
989:       ETOTAL=USEFRAC*ETOTALTMP+(1.0D0-USEFRAC)*ETOTAL976:       ETOTAL=USEFRAC*ETOTALTMP+(1.0D0-USEFRAC)*ETOTAL
990:       G(1:D)=USEFRAC*MYGTMP(1:D)+(1.0D0-USEFRAC)*G(1:D)977:       G(1:D)=USEFRAC*MYGTMP(1:D)+(1.0D0-USEFRAC)*G(1:D)
1001:       WRITE(MYUNIT,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/INTIMAGE,COLDFUSIONLIMIT988:       WRITE(MYUNIT,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/INTIMAGE,COLDFUSIONLIMIT
1002:       DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT)989:       DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT)
1003:       DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &990:       DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &
1004:   &              DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)991:   &              DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)
1005:       INTIMAGE=INTIMAGESAVE992:       INTIMAGE=INTIMAGESAVE
1006:       RETURN993:       RETURN
1007:    ENDIF994:    ENDIF
1008: 995: 
1009:    STEPTOT = SUM(STEPIMAGE)/INTIMAGE996:    STEPTOT = SUM(STEPIMAGE)/INTIMAGE
1010: 997: 
1011:    MAXRMS=-1.0D0 
1012:    MAXEEE=-1.0D100 
1013:    DO J1=2,INTIMAGE+1 
1014:       IF (EEE(J1).GT.MAXEEE) THEN 
1015:          MAXEEE=EEE(J1) 
1016:          JMAXEEE=J1 
1017:       ENDIF 
1018:       DUMMY=0.0D0 
1019:       DO J2=1,3*NATOMS 
1020:          DUMMY=DUMMY+GGG(3*NATOMS*(J1-1)+J2)**2 
1021:       ENDDO 
1022:       IF (DUMMY.GT.MAXRMS) THEN 
1023:          MAXRMS=DUMMY 
1024:          JMAXRMS=J1 
1025:       ENDIF 
1026:    ENDDO 
1027:    MAXRMS=SQRT(MAXRMS/(3*NACTIVE)) 
1028:  
1029:    IF (DEBUG) THEN998:    IF (DEBUG) THEN
1030:       WRITE(MYUNIT,'(A,I6,2G20.10,3(G20.10,I8))') ' intlbfgs> steps: ',NITERDONE,ETOTAL/INTIMAGE,RMS,STEPTOT,NACTIVE, &999:       WRITE(MYUNIT,'(A,I6,2G20.10,G20.10,I8)') ' intlbfgs> steps: ',NITERDONE,ETOTAL/INTIMAGE,RMS,STEPTOT,NACTIVE
1031:   &                                                        MAXEEE,JMAXEEE,MAXRMS,JMAXRMS 
1032:       CALL FLUSH(6)1000:       CALL FLUSH(6)
1033:    ENDIF1001:    ENDIF
1034: 1002: 
1035:    IF (.NOT.SWITCHED) THEN1003:    IF (.NOT.SWITCHED) THEN
1036: !     IF ((NITERDONE-NLASTGOODE.GT.INTRELSTEPS).AND.((ETOTAL.GT.LASTGOODE).OR.(ETOTAL/INTIMAGE.GT.MAXCONE*1.0D8))) THEN1004: !     IF ((NITERDONE-NLASTGOODE.GT.INTRELSTEPS).AND.((ETOTAL.GT.LASTGOODE).OR.(ETOTAL/INTIMAGE.GT.MAXCONE*1.0D8))) THEN
1037:       IF (.FALSE.) THEN ! no backtracking1005:       IF (.FALSE.) THEN ! no backtracking
1038:          WRITE(MYUNIT,'(2(A,I6))') ' intlbfgs> Backtracking ',NBACKTRACK,' steps, current active atoms=',NACTIVE1006:          WRITE(MYUNIT,'(2(A,I6))') ' intlbfgs> Backtracking ',NBACKTRACK,' steps, current active atoms=',NACTIVE
1039:          NTRIES(NEWATOM)=NTRIES(NEWATOM)+11007:          NTRIES(NEWATOM)=NTRIES(NEWATOM)+1
1040:          IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.1008:          IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.
1041: !1009: !
1122:          IF (NDUMMY.NE.NACTIVE) THEN1090:          IF (NDUMMY.NE.NACTIVE) THEN
1123:             WRITE(MYUNIT,'(A,I6)') ' intlbfgs> ERROR *** inconsistency in number of active atoms. ',NDUMMY,' should be ',NACTIVE1091:             WRITE(MYUNIT,'(A,I6)') ' intlbfgs> ERROR *** inconsistency in number of active atoms. ',NDUMMY,' should be ',NACTIVE
1124:             DO J1=1,NATOMS1092:             DO J1=1,NATOMS
1125:                IF (ATOMACTIVE(J1)) WRITE(MYUNIT,'(A,I6)') ' active atom ',J11093:                IF (ATOMACTIVE(J1)) WRITE(MYUNIT,'(A,I6)') ' active atom ',J1
1126:             ENDDO1094:             ENDDO
1127:             STOP1095:             STOP
1128:          ENDIF1096:          ENDIF
1129:          ADDATOM=.TRUE.1097:          ADDATOM=.TRUE.
1130: 1098: 
1131:          CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)1099:          CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
1132:          IF (QCIADDREP.GT.0) THEN1100:          IF (CHECKCONINT) THEN
1133:             CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
1134:          ELSEIF (CHECKCONINT) THEN 
1135:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1101:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1136:          ELSE1102:          ELSE
1137:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1103:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1138:          ENDIF1104:          ENDIF
1139:       ENDIF1105:       ENDIF
1140:       LASTGOODE=ETOTAL1106:       LASTGOODE=ETOTAL
1141:    ENDIF1107:    ENDIF
1142:  
1143:    EXITSTATUS=01108:    EXITSTATUS=0
1144:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW1109:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW
1145: !  IF ((.NOT.SWITCHED).AND.(RMS<=INTRMSTOL*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))).AND.NITERDONE>1) EXITSTATUS=1 1110:    IF ((.NOT.SWITCHED).AND.(RMS<=INTRMSTOL*NACTIVE/NATOMS).AND.NITERDONE>1) EXITSTATUS=1 
1146:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 1111:    IF (SWITCHED.AND.(RMS<=CQMAX).AND.NITERDONE>1) EXITSTATUS=1 
1147: !  IF (SWITCHED.AND.(RMS<=CQMAX).AND.NITERDONE>1) EXITSTATUS=1  
1148:    IF (SWITCHED.AND.(MAXRMS<=CQMAX).AND.NITERDONE>1) EXITSTATUS=1  
1149:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=21112:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=2
1150: !  IF (SQRT(ADMAX).GT.IMSEPMAX) EXITSTATUS=0 ! prevent converge if largest atomic displacement is too big1113: !  IF (SQRT(ADMAX).GT.IMSEPMAX) EXITSTATUS=0 ! prevent converge if largest atomic displacement is too big
1151:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW1114:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW
1152: 1115: 
1153:    IF (EXITSTATUS > 0) THEN  1116:    IF (EXITSTATUS > 0) THEN  
1154:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on1117:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on
1155: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771118:          IF (ETOTAL/INTIMAGE.GT.MAXCONE*NACTIVE/NATOMS) GOTO 777
1156:          IF (MAXEEE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777 
1157:          IF (NACTIVE.LT.NATOMS) THEN 1119:          IF (NACTIVE.LT.NATOMS) THEN 
1158:             ADDATOM=.TRUE.1120:             ADDATOM=.TRUE.
1159:             GOTO 7771121:             GOTO 777
1160:          ENDIF1122:          ENDIF
1161:          CALL MYCPU_TIME(FTIME,.FALSE.)1123:          CALL MYCPU_TIME(FTIME,.FALSE.)
1162:          WRITE(MYUNIT,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &1124:          WRITE(MYUNIT,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &
1163:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME1125:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME
1164:          IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ)1126:          IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ)
1165:          IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE,MYUNIT)1127:          IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE,MYUNIT)
1166:          WRITE(MYUNIT,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'1128:          WRITE(MYUNIT,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'
1319:   &      DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)1281:   &      DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)
1320: INTIMAGE=INTIMAGESAVE1282: INTIMAGE=INTIMAGESAVE
1321: 1283: 
1322: STOP1284: STOP
1323: 1285: 
1324: END SUBROUTINE INTLBFGS1286: END SUBROUTINE INTLBFGS
1325: !1287: !
1326: ! Neighbour list for repulsions to reduce cost of constraint potential.1288: ! Neighbour list for repulsions to reduce cost of constraint potential.
1327: !1289: !
1328: SUBROUTINE CHECKREP(INTIMAGE,XYZ,NOPT,NNSTART,NSTART)1290: SUBROUTINE CHECKREP(INTIMAGE,XYZ,NOPT,NNSTART,NSTART)
1329: USE COMMONS,ONLY : NREPI, NREPJ, NREPCUT, NNREPULSIVE, NREPULSIVE, REPI, REPJ, REPCUT, CHECKREPCUTOFF, DEBUG, MYUNIT, &1291: USE COMMONS,ONLY : NREPI, NREPJ, NREPCUT, NNREPULSIVE, NREPULSIVE, REPI, REPJ, REPCUT, CHECKREPCUTOFF, DEBUG, MYUNIT, INTFROZEN
1330:   &                INTFROZEN, NNREPULSIVE, intconstraintrep 
1331: USE PORFUNCS1292: USE PORFUNCS
1332: IMPLICIT NONE1293: IMPLICIT NONE
1333: INTEGER JJ, KK, NI1, NJ1, NI2, NJ2, INTIMAGE, NOPT, NI, NJ, NNSTART, NSTART1294: INTEGER JJ, KK, NI1, NJ1, NI2, NJ2, INTIMAGE, NOPT, NI, NJ, NNSTART, NSTART
1334: DOUBLE PRECISION LDIST, XYZ(NOPT*(INTIMAGE+2)),COMPARE1295: DOUBLE PRECISION LDIST, XYZ(NOPT*(INTIMAGE+2)),COMPARE
1335: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,DMIN1296: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,DMIN
1336: LOGICAL NOINT1297: LOGICAL NOINT
1337: 1298: 
1338: IF (INTCONSTRAINTREP.EQ.0) THEN 
1339:    NNREPULSIVE=0 
1340:    RETURN 
1341: ENDIF 
1342:  
1343: NNREPULSIVE=NNSTART1299: NNREPULSIVE=NNSTART
1344: DO JJ=NSTART,NREPULSIVE1300: DO JJ=NSTART,NREPULSIVE
1345:    COMPARE=(CHECKREPCUTOFF*REPCUT(JJ))**21301:    COMPARE=(CHECKREPCUTOFF*REPCUT(JJ))**2
1346:    NI=REPI(JJ)1302:    NI=REPI(JJ)
1347:    NJ=REPJ(JJ)1303:    NJ=REPJ(JJ)
1348:    DO KK=1,INTIMAGE+2 ! first check for standard distances within threshold1304:    DO KK=1,INTIMAGE+2 ! first check for standard distances within threshold
1349:       LDIST=(XYZ((KK-1)*NOPT+3*(NI-1)+1)-XYZ((KK-1)*NOPT+3*(NJ-1)+1))**2 &1305:       LDIST=(XYZ((KK-1)*NOPT+3*(NI-1)+1)-XYZ((KK-1)*NOPT+3*(NJ-1)+1))**2 &
1350:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+2)-XYZ((KK-1)*NOPT+3*(NJ-1)+2))**2 &1306:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+2)-XYZ((KK-1)*NOPT+3*(NJ-1)+2))**2 &
1351:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+3)-XYZ((KK-1)*NOPT+3*(NJ-1)+3))**21307:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+3)-XYZ((KK-1)*NOPT+3*(NJ-1)+3))**2
1352:       IF (LDIST.LT.COMPARE) THEN1308:       IF (LDIST.LT.COMPARE) THEN
1446: CLOSE(UNIT)1402: CLOSE(UNIT)
1447: WRITE(MYUNIT,'(A)') ' writeprofile> Interpolated energy profile was saved to file "'//trim(filename)//'"'1403: WRITE(MYUNIT,'(A)') ' writeprofile> Interpolated energy profile was saved to file "'//trim(filename)//'"'
1448: 1404: 
1449: END SUBROUTINE WRITEPROFILE1405: END SUBROUTINE WRITEPROFILE
1450: 1406: 
1451: SUBROUTINE DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE)1407: SUBROUTINE DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE)
1452: USE COMMONS, ONLY : CONACTIVE, CONI, CONJ, ATOMACTIVE, CONDISTREF, REPI, REPJ, REPCUT, INTREPSEP,  &1408: USE COMMONS, ONLY : CONACTIVE, CONI, CONJ, ATOMACTIVE, CONDISTREF, REPI, REPJ, REPCUT, INTREPSEP,  &
1453:   &             INTCONSTRAINREPCUT, NREPULSIVE, NREPMAX, MAXCONUSE, CHECKCONINT, &1409:   &             INTCONSTRAINREPCUT, NREPULSIVE, NREPMAX, MAXCONUSE, CHECKCONINT, &
1454:   &             FREEZENODEST, NNREPULSIVE, NATOMS, DEBUG, MYUNIT, INTFROZEN, &1410:   &             FREEZENODEST, NNREPULSIVE, NATOMS, DEBUG, MYUNIT, INTFROZEN, &
1455:   &             NREPULSIVEFIX, REPIFIX, REPJFIX, REPCUTFIX, NREPI, NREPJ, NREPCUT, MAXNACTIVE, &1411:   &             NREPULSIVEFIX, REPIFIX, REPJFIX, REPCUTFIX, NREPI, NREPJ, NREPCUT, MAXNACTIVE, &
1456:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, QCIADDREP1412:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT
1457: IMPLICIT NONE1413: IMPLICIT NONE
1458: INTEGER INTIMAGE1414: INTEGER INTIMAGE
1459: INTEGER NBEST, NCONTOACTIVE(NATOMS),  NCONSTRAINT, J2, NTRIES(NATOMS), NEWATOM,  CONLIST(NATOMS), N1, N2, N3, &1415: INTEGER NBEST, NCONTOACTIVE(NATOMS),  NCONSTRAINT, J2, NTRIES(NATOMS), NEWATOM,  CONLIST(NATOMS), N1, N2, N3, &
1460:   &     NTOADD, NADDED, NMININT, NMAXINT, TURNONORDER(NATOMS), NDUMMY, J1, J3, NITERDONE, NCONFORNEWATOM, NACTIVE1416:   &     NTOADD, NADDED, NMININT, NMAXINT, TURNONORDER(NATOMS), NDUMMY, J1, J3, NITERDONE, NCONFORNEWATOM, NACTIVE
1461: DOUBLE PRECISION DUMMY, DUMMY2, DPRAND, RANDOM, CONDIST(NATOMS), DMIN1417: DOUBLE PRECISION DUMMY, DUMMY2, DPRAND, RANDOM, CONDIST(NATOMS), DMIN
1462: INTEGER NDFORNEWATOM, BESTPRESERVEDN(NATOMS)1418: INTEGER NDFORNEWATOM, BESTPRESERVEDN(NATOMS)
1463: DOUBLE PRECISION BESTPRESERVEDD(NATOMS), BESTCLOSESTD(NATOMS), INVDTOACTIVE(NATOMS)1419: DOUBLE PRECISION BESTPRESERVEDD(NATOMS), BESTCLOSESTD(NATOMS), INVDTOACTIVE(NATOMS)
1464: LOGICAL IMGFREEZE(INTIMAGE), ADDREP(NATOMS)1420: LOGICAL IMGFREEZE(INTIMAGE), ADDREP(NATOMS)
1465: DOUBLE PRECISION C1, C2, C3, VEC1(3), VEC2(3), VEC3(3), ESAVED, ESAVEC, ESAVE01421: DOUBLE PRECISION C1, C2, C3, VEC1(3), VEC2(3), VEC3(3), ESAVED, ESAVEC, ESAVE0
1466: INTEGER NCFORNEWATOM, BESTCLOSESTN(NATOMS), NNREPSAVE, NREPSAVE1422: INTEGER NCFORNEWATOM, BESTCLOSESTN(NATOMS), NNREPSAVE, NREPSAVE
1561:             IF (.NOT.ATOMACTIVE(CONJ(J1))) THEN1517:             IF (.NOT.ATOMACTIVE(CONJ(J1))) THEN
1562:                IF (CONDISTREF(J1).LT.DUMMY2) THEN1518:                IF (CONDISTREF(J1).LT.DUMMY2) THEN
1563:                   DUMMY2=CONDISTREF(J1)1519:                   DUMMY2=CONDISTREF(J1)
1564:                   NEWATOM=CONJ(J1)1520:                   NEWATOM=CONJ(J1)
1565:                ENDIF1521:                ENDIF
1566:             ENDIF1522:             ENDIF
1567:          ENDIF1523:          ENDIF
1568:       ENDDO1524:       ENDDO
1569:       IF (DEBUG) WRITE(MYUNIT,'(3(A,I6),A,F15.5)') ' intlbfgs> Choosing new active atom ',NEWATOM,' new constraints=', &1525:       IF (DEBUG) WRITE(MYUNIT,'(3(A,I6),A,F15.5)') ' intlbfgs> Choosing new active atom ',NEWATOM,' new constraints=', &
1570:   &                                       NCONTOACTIVE(NEWATOM),' maximum=',NBEST,' shortest constraint=',DUMMY21526:   &                                       NCONTOACTIVE(NEWATOM),' maximum=',NBEST,' shortest constraint=',DUMMY2
1571: !     IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ) 
1572: !     IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE,MYUNIT) 
1573:  
1574:           1527:           
1575:       IF (NEWATOM*NBEST.EQ.0) THEN ! sanity check1528:       IF (NEWATOM*NBEST.EQ.0) THEN ! sanity check
1576:          WRITE(MYUNIT,'(A,I6,A,2I6)') ' intlbfgs> ERROR *** new active atom not set'1529:          WRITE(MYUNIT,'(A,I6,A,2I6)') ' intlbfgs> ERROR *** new active atom not set'
1577:          STOP1530:          STOP
1578:       ELSE1531:       ELSE
1579: !1532: !
1580: !  We need a sorted list of up to 3 active atoms, sorted according to how well the1533: !  We need a sorted list of up to 3 active atoms, sorted according to how well the
1581: !  end point distance is preserved, even if they don't satisfy the constraint 1534: !  end point distance is preserved, even if they don't satisfy the constraint 
1582: !  condition. We want three atoms to use for a local axis system in the interpolation.1535: !  condition. We want three atoms to use for a local axis system in the interpolation.
1583: !1536: !
1878:                VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)1831:                VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)
1879:                DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)1832:                DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)
1880:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY1833:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY
1881:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)1834:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
1882:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)1835:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)
1883:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)1836:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
1884:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &1837:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &
1885:   &            XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)+C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)1838:   &            XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)+C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)
1886:             ENDDO1839:             ENDDO
1887:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list1840:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
1888:             IF (QCIADDREP.GT.0) THEN1841:             IF (CHECKCONINT) THEN
1889:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
1890:             ELSEIF (CHECKCONINT) THEN 
1891:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1842:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1892:             ELSE1843:             ELSE
1893:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1844:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1894:             ENDIF1845:             ENDIF
1895:             ESAVE0=ETOTAL1846:             ESAVE0=ETOTAL
1896:             DO J1=2,INTIMAGE+11847:             DO J1=2,INTIMAGE+1
1897:                XSAVE0(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)1848:                XSAVE0(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
1898:             ENDDO1849:             ENDDO
1899:          ENDIF1850:          ENDIF
1900:          IF (NDFORNEWATOM.GE.3) THEN1851:          IF (NDFORNEWATOM.GE.3) THEN
1983:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &1934:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &
1984:   &            XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3)+ &1935:   &            XYZ((J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTPRESERVEDN(N1)-1)+3)+ &
1985:   &                   C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)+0.01D0*(DPRAND()-0.5D0)*2.0D01936:   &                   C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)+0.01D0*(DPRAND()-0.5D0)*2.0D0
1986: !              WRITE(MYUNIT,'(A,I6,3G20.10)') 'intlbfgs> J1,C1,C2,C3=',J1,C1,C2,C31937: !              WRITE(MYUNIT,'(A,I6,3G20.10)') 'intlbfgs> J1,C1,C2,C3=',J1,C1,C2,C3
1987: !              WRITE(MYUNIT,'(A,9G20.10)') 'intlbfgs> VEC1,2,3=',VEC1(1:3),VEC2(1:3),VEC3(1:3)1938: !              WRITE(MYUNIT,'(A,9G20.10)') 'intlbfgs> VEC1,2,3=',VEC1(1:3),VEC2(1:3),VEC3(1:3)
1988: !              WRITE(MYUNIT,'(A,6I6)') 'intlbfgs> N1,N2,N3,Bestpreserved N1,N2,N3=',N1,N2,N3, &1939: !              WRITE(MYUNIT,'(A,6I6)') 'intlbfgs> N1,N2,N3,Bestpreserved N1,N2,N3=',N1,N2,N3, &
1989: ! &                 BESTPRESERVEDN(N1),BESTPRESERVEDN(N2),BESTPRESERVEDN(N3)1940: ! &                 BESTPRESERVEDN(N1),BESTPRESERVEDN(N2),BESTPRESERVEDN(N3)
1990:             ENDDO1941:             ENDDO
1991: 1942: 
1992:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list1943:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
1993:             IF (QCIADDREP.GT.0) THEN1944:             IF (CHECKCONINT) THEN
1994:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
1995:             ELSEIF (CHECKCONINT) THEN 
1996:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1945:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1997:             ELSE1946:             ELSE
1998:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1947:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1999:             ENDIF1948:             ENDIF
2000:             ESAVED=ETOTAL1949:             ESAVED=ETOTAL
2001:             DO J1=2,INTIMAGE+11950:             DO J1=2,INTIMAGE+1
2002:                XSAVED(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)1951:                XSAVED(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2003:             ENDDO1952:             ENDDO
2004:          ENDIF1953:          ENDIF
2005: 1954: 
2073:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY2022:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY
2074:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)2023:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
2075:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)2024:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)
2076:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)2025:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
2077:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &2026:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &
2078:   &            XYZ((J1-1)*3*NATOMS+3*(BESTCLOSESTN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTCLOSESTN(N1)-1)+3)+ &2027:   &            XYZ((J1-1)*3*NATOMS+3*(BESTCLOSESTN(N1)-1)+1:(J1-1)*3*NATOMS+3*(BESTCLOSESTN(N1)-1)+3)+ &
2079:   &                   C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)2028:   &                   C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)
2080:             ENDDO2029:             ENDDO
2081: 2030: 
2082:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2031:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2083:             IF (QCIADDREP.GT.0) THEN2032:             IF (CHECKCONINT) THEN
2084:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2085:             ELSEIF (CHECKCONINT) THEN 
2086:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2033:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2087:             ELSE2034:             ELSE
2088:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2035:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2089:             ENDIF2036:             ENDIF
2090:             ESAVEC=ETOTAL2037:             ESAVEC=ETOTAL
2091:             DO J1=2,INTIMAGE+12038:             DO J1=2,INTIMAGE+1
2092:                XSAVEC(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2039:                XSAVEC(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2093:             ENDDO2040:             ENDDO
2094:          ENDIF2041:          ENDIF
2095: !2042: !
2104:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))/(INTIMAGE+1) &2051:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))/(INTIMAGE+1) &
2105:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+1)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+1))/(INTIMAGE+1)2052:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+1)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+1))/(INTIMAGE+1)
2106:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+2)  &2053:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+2)  &
2107:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))/(INTIMAGE+1) &2054:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))/(INTIMAGE+1) &
2108:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+2)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+2))/(INTIMAGE+1)2055:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+2)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+2))/(INTIMAGE+1)
2109:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)  &2056:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)  &
2110:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))/(INTIMAGE+1) &2057:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))/(INTIMAGE+1) &
2111:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+3)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+3))/(INTIMAGE+1)2058:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+3)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+3))/(INTIMAGE+1)
2112:          ENDDO2059:          ENDDO
2113:          CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2060:          CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2114:          IF (QCIADDREP.GT.0) THEN2061:          IF (CHECKCONINT) THEN
2115:             CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2116:          ELSEIF (CHECKCONINT) THEN 
2117:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2062:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2118:          ELSE2063:          ELSE
2119:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2064:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2120:          ENDIF2065:          ENDIF
2121:          IF (DEBUG) WRITE(MYUNIT,'(A,4G15.5)') ' intlbfgs> energies for constrained, preserved, closest, and linear schemes=', &2066:          IF (DEBUG) WRITE(MYUNIT,'(A,4G15.5)') ' intlbfgs> energies for constrained, preserved, closest, and linear schemes=', &
2122:   &                 ESAVE0,ESAVED,ESAVEC,ETOTAL2067:   &                 ESAVE0,ESAVED,ESAVEC,ETOTAL
2123:          IF ((ETOTAL.LT.ESAVEC).AND.(ETOTAL.LT.ESAVED).AND.(ETOTAL.LT.ESAVE0)) THEN2068:          IF ((ETOTAL.LT.ESAVEC).AND.(ETOTAL.LT.ESAVED).AND.(ETOTAL.LT.ESAVE0)) THEN
2124:             IF (DEBUG) WRITE(MYUNIT,'(A,2G20.10)') ' intlbfgs> lowest energy from linear interpolation'2069:             IF (DEBUG) WRITE(MYUNIT,'(A,2G20.10)') ' intlbfgs> lowest energy from linear interpolation'
2125:          ELSE IF ((ESAVEC.LT.ESAVED).AND.(ESAVEC.LT.ESAVE0)) THEN2070:          ELSE IF ((ESAVEC.LT.ESAVED).AND.(ESAVEC.LT.ESAVE0)) THEN
2126:             IF (DEBUG) WRITE(MYUNIT,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest atoms'2071:             IF (DEBUG) WRITE(MYUNIT,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest atoms'
2172:       ENDIF2117:       ENDIF
2173: !2118: !
2174: ! Turn frozen images off for new added atom.2119: ! Turn frozen images off for new added atom.
2175: !2120: !
2176: !     IF (DEBUG) WRITE(MYUNIT,'(A)') ' intlbfgs> turning off frozen images'2121: !     IF (DEBUG) WRITE(MYUNIT,'(A)') ' intlbfgs> turning off frozen images'
2177: !     IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.2122: !     IF (FREEZENODEST) IMGFREEZE(1:INTIMAGE)=.FALSE.
2178:       CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2123:       CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2179: !2124: !
2180: ! need a new gradient since the active atom has changed !2125: ! need a new gradient since the active atom has changed !
2181: !2126: !
2182:       IF (QCIADDREP.GT.0) THEN2127:       IF (CHECKCONINT) THEN
2183:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2184:       ELSEIF (CHECKCONINT) THEN 
2185:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2128:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2186:       ELSE2129:       ELSE
2187:          CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2130:          CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2188:       ENDIF2131:       ENDIF
2189: 2132: 
2190: END SUBROUTINE DOADDATOM2133: END SUBROUTINE DOADDATOM
2191: 2134: 
2192: SUBROUTINE CHECKPERC(LXYZ,LINTCONSTRAINTTOL,NQCIFREEZE,NCPFIT)2135: SUBROUTINE CHECKPERC(LXYZ,LINTCONSTRAINTTOL,NQCIFREEZE,NCPFIT)
2193: USE COMMONS, ONLY : ATOMACTIVE, NCONSTRAINT, INTFROZEN, CONI, CONJ, CONDISTREF, INTCONMAX, INTCONSTRAINTTOL, &2136: USE COMMONS, ONLY : ATOMACTIVE, NCONSTRAINT, INTFROZEN, CONI, CONJ, CONDISTREF, INTCONMAX, INTCONSTRAINTTOL, &
2194:   &             INTCONSEP, NCONGEOM, CONGEOM, CONIFIX, CONJFIX, CONDISTREFFIX, MYUNIT, INTCONCUT, &2137:   &             INTCONSEP, NCONGEOM, CONGEOM, CONIFIX, CONJFIX, CONDISTREFFIX, MYUNIT, INTCONCUT, &
2195:   &             NCONSTRAINTFIX, PERIODIC, TWOD, RIGID, CONDATT, CONCUT, CONCUTFIX, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, &2138:   &             NCONSTRAINTFIX, PERIODIC, TWOD, RIGID, CONDATT, CONCUT, CONCUTFIX, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, &
2196:   &             BONDS, QCIAMBERT, QCIADDREP, QCIADDREPCUT, QCIBONDS, QCISECOND2139:   &             BONDS, QCIAMBERT
2197: IMPLICIT NONE2140: IMPLICIT NONE
2198: INTEGER NDIST1(NATOMS), NCYCLE, DMIN1, DMAX1, NUNCON1, J1, J2, J3, NQCIFREEZE, J4, NCPFIT, LUNIT, GETUNIT2141: INTEGER NDIST1(NATOMS), NCYCLE, DMIN1, DMAX1, NUNCON1, J1, J2, J3, NQCIFREEZE, J4, NCPFIT, LUNIT, GETUNIT
2199: INTEGER NI1, NJ1, NI2, NJ2, J5, ATOM1, ATOM2 
2200: DOUBLE PRECISION LINTCONSTRAINTTOL, MAXCONDIST, MINCONDIST, DS, DF, LXYZ((3*NATOMS)*2)2142: DOUBLE PRECISION LINTCONSTRAINTTOL, MAXCONDIST, MINCONDIST, DS, DF, LXYZ((3*NATOMS)*2)
2201: DOUBLE PRECISION DSMIN, DSMAX, DSMEAN, D, DIST2, RMAT(3,3), DUMMY, X1, Y1, Z1, X2, Y2, Z2, DMIN, D22143: DOUBLE PRECISION DSMIN, DSMAX, DSMEAN, D, DIST2, RMAT(3,3), DUMMY
2202: LOGICAL CHANGED, LDEBUG, CONFILET2144: LOGICAL CHANGED, LDEBUG, CONFILET
2203: LOGICAL :: CALLED=.FALSE.2145: LOGICAL :: CALLED=.FALSE.
2204: SAVE CALLED2146: SAVE CALLED
2205: !for QCIAMBER2147: !for QCIAMBER
2206: INTEGER NBOND2148: INTEGER NBOND
2207: 2149: 
2208: LINTCONSTRAINTTOL=INTCONSTRAINTTOL2150: LINTCONSTRAINTTOL=INTCONSTRAINTTOL
2209: 2151: 
2210: IF (.NOT.ALLOCATED(ATOMACTIVE)) ALLOCATE(ATOMACTIVE(NATOMS))2152: IF (.NOT.ALLOCATED(ATOMACTIVE)) ALLOCATE(ATOMACTIVE(NATOMS))
2211: !2153: !
2260:       ENDDO2202:       ENDDO
2261:    ENDIF2203:    ENDIF
2262:    ALLOCATE(CONIFIX(INTCONMAX),CONJFIX(INTCONMAX),CONCUTFIX(INTCONMAX),CONDISTREFFIX(INTCONMAX))2204:    ALLOCATE(CONIFIX(INTCONMAX),CONJFIX(INTCONMAX),CONCUTFIX(INTCONMAX),CONDISTREFFIX(INTCONMAX))
2263: ENDIF2205: ENDIF
2264: 2206: 
2265: INQUIRE(FILE='constraintfile',EXIST=CONFILET)2207: INQUIRE(FILE='constraintfile',EXIST=CONFILET)
2266: 2208: 
2267: 51   NCONSTRAINT=0 2209: 51   NCONSTRAINT=0 
2268: MAXCONDIST=-1.0D02210: MAXCONDIST=-1.0D0
2269: MINCONDIST=1.0D1002211: MINCONDIST=1.0D100
2270: IF (QCIAMBERT) THEN             2212: IF (QCIAMBERT) THEN             !kr366> assume we use two endpoints and topology for amber constraints
2271:    CALL TOPOLOGY_READER(NBOND)  2213:    CALL TOPOLOGY_READER(NBOND)  !get number of bonds and bonds from topology
2272: ! 
2273: !  kr366> assume we use two endpoints and topology for amber constraints 
2274: !  get number of bonds and bonds from topology 
2275: !  loop through all bonds and add them to constraint list 
2276: ! 
2277:    DO J2=1,NBOND                !loop through all bonds and add them to constraint list2214:    DO J2=1,NBOND                !loop through all bonds and add them to constraint list
2278:       IF (INTFROZEN(BONDS(J2,1)).AND.INTFROZEN(BONDS(J2,2))) CYCLE ! no constraints between intfrozen atoms2215:       IF (INTFROZEN(BONDS(J2,1)).AND.INTFROZEN(BONDS(J2,2))) THEN
2279:       NCONSTRAINT=NCONSTRAINT+12216:          INTFROZEN(BONDS(J2,1))=.FALSE.
2280:       WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for atoms ',BONDS(J2,1),BONDS(J2,2), &2217:          INTFROZEN(BONDS(J2,2))=.FALSE.
2281:   &                     '  total=',NCONSTRAINT2218:          WRITE(MYUNIT,'(A,2I6)') 'qciamber> Unfreeze bonded atoms: ', BONDS(J2,1),BONDS(J2,2)
 2219:          NQCIFREEZE=NQCIFREEZE-2
 2220:       ENDIF
2282:       DS=SQRT((LXYZ(3*(BONDS(J2,1)-1)+1)-LXYZ(3*(BONDS(J2,2)-1)+1))**2 &2221:       DS=SQRT((LXYZ(3*(BONDS(J2,1)-1)+1)-LXYZ(3*(BONDS(J2,2)-1)+1))**2 &
2283:   &          +(LXYZ(3*(BONDS(J2,1)-1)+2)-LXYZ(3*(BONDS(J2,2)-1)+2))**2 &2222:   &          +(LXYZ(3*(BONDS(J2,1)-1)+2)-LXYZ(3*(BONDS(J2,2)-1)+2))**2 &
2284:   &          +(LXYZ(3*(BONDS(J2,1)-1)+3)-LXYZ(3*(BONDS(J2,2)-1)+3))**2)2223:   &          +(LXYZ(3*(BONDS(J2,1)-1)+3)-LXYZ(3*(BONDS(J2,2)-1)+3))**2)
2285:       DF=SQRT((LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+1)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+1))**2 &2224:       DF=SQRT((LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+1)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+1))**2 &
2286:   &          +(LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+2)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+2))**2 &2225:   &          +(LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+2)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+2))**2 &
2287:   &          +(LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+3)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+3))**2)2226:   &          +(LXYZ(3*NATOMS+3*(BONDS(J2,1)-1)+3)-LXYZ(3*NATOMS+3*(BONDS(J2,2)-1)+3))**2)
2288:       IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE2227:       WRITE(MYUNIT,'(A,2I6,2G20.10)') 'intlbfgs> J2,J3,DS,DF=', BONDS(J2,1),BONDS(J2,2),DS,DF
2289:       CONI(NCONSTRAINT)=MIN(BONDS(J2,1),BONDS(J2,2))2228:       NCONSTRAINT=NCONSTRAINT+1
2290:       CONJ(NCONSTRAINT)=MAX(BONDS(J2,2),BONDS(J2,2))2229:       WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for atoms ',BONDS(J2,1),BONDS(J2,2), &
 2230:   &                     '  total=',NCONSTRAINT
 2231:       CONI(NCONSTRAINT)=BONDS(J2,1)
 2232:       CONJ(NCONSTRAINT)=BONDS(J2,2)
2291:       CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D02233:       CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D0
2292:       CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D02234:       CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D0
2293:       IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)2235:       IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)
2294:       IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)2236:       IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)
2295:    ENDDO2237:    ENDDO
2296:    WRITE(MYUNIT,'(A,2I6,A,3F12.2,A,F12.4,A,I8)') ' intlbfgs> constrain distance for atoms ',CONI(NCONSTRAINT), &2238:    IF (DEBUG) WRITE(MYUNIT,'(A,I10,2(A,G20.10))') ' intlbfgs> Total distance constraints=',NCONSTRAINT, &
2297:   &             CONJ(NCONSTRAINT),' values are ',DS,DF,DUMMY,' fraction=',2*ABS(DS-DF)/(DS+DF), &2239:   &                                     ' shortest=',MINCONDIST,' longest=',MAXCONDIST      
2298:   &            ' # bond constraints=',NCONSTRAINT2240:   
2299:    QCIBONDS=NCONSTRAINT2241: ELSE IF (CONFILET) THEN 
2300: ! 
2301: ! Add constraints for second-nearest neighbours - should correspond to bond angles 
2302: ! 
2303:    DO J2=1,NBOND 
2304:       inloop: DO J3=J2+1,NBOND 
2305:         IF (BONDS(J2,1).EQ.BONDS(J3,1)) THEN 
2306:            ATOM1=BONDS(J2,2) 
2307:            ATOM2=BONDS(J3,2) 
2308:         ELSEIF (BONDS(J2,1).EQ.BONDS(J3,2)) THEN 
2309:            ATOM1=BONDS(J2,2) 
2310:            ATOM2=BONDS(J3,1) 
2311:         ELSEIF (BONDS(J2,2).EQ.BONDS(J3,1)) THEN 
2312:            ATOM1=BONDS(J2,1) 
2313:            ATOM2=BONDS(J3,2) 
2314:         ELSEIF (BONDS(J2,2).EQ.BONDS(J3,2)) THEN 
2315:            ATOM1=BONDS(J2,1) 
2316:            ATOM2=BONDS(J3,1) 
2317:         ELSE 
2318:            CYCLE inloop 
2319:         ENDIF 
2320:         IF (INTFROZEN(ATOM1).AND.INTFROZEN(ATOM2)) CYCLE ! no constraints between intfrozen atoms 
2321:         NCONSTRAINT=NCONSTRAINT+1 
2322:         WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for second neighbours ',ATOM1,ATOM2, & 
2323:   &                     '  total=',NCONSTRAINT 
2324:          DS=SQRT((LXYZ(3*(ATOM1-1)+1)-LXYZ(3*(ATOM2-1)+1))**2 & 
2325:   &             +(LXYZ(3*(ATOM1-1)+2)-LXYZ(3*(ATOM2-1)+2))**2 & 
2326:   &             +(LXYZ(3*(ATOM1-1)+3)-LXYZ(3*(ATOM2-1)+3))**2) 
2327:          DF=SQRT((LXYZ(3*NATOMS+3*(ATOM1-1)+1)-LXYZ(3*NATOMS+3*(ATOM2-1)+1))**2 & 
2328:   &             +(LXYZ(3*NATOMS+3*(ATOM1-1)+2)-LXYZ(3*NATOMS+3*(ATOM2-1)+2))**2 & 
2329:   &             +(LXYZ(3*NATOMS+3*(ATOM1-1)+3)-LXYZ(3*NATOMS+3*(ATOM2-1)+3))**2) 
2330:          IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE 
2331:          CONI(NCONSTRAINT)=MIN(ATOM1,ATOM2) 
2332:          CONJ(NCONSTRAINT)=MAX(ATOM1,ATOM2) 
2333:          CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D0 
2334:          CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D0 
2335:          IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT) 
2336:          IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT) 
2337:          WRITE(MYUNIT,'(A,2I6,A,3F12.2,A,F12.4,A,2I8)') ' intlbfgs> constrain distance for atoms ',CONI(NCONSTRAINT), & 
2338:   &             CONJ(NCONSTRAINT),' values are ',DS,DF,DUMMY,' fraction=',2*ABS(DS-DF)/(DS+DF), & 
2339:   &            ' # second neighbour constraints and total=',QCISECOND,NCONSTRAINT 
2340:       ENDDO inloop 
2341:    ENDDO 
2342:    QCISECOND=NCONSTRAINT-QCIBONDS 
2343: ENDIF 
2344: IF (CONFILET) THEN  
2345:     LUNIT=GETUNIT()2242:     LUNIT=GETUNIT()
2346:     OPEN(LUNIT,FILE='constraintfile',STATUS='OLD')2243:     OPEN(LUNIT,FILE='constraintfile',STATUS='OLD')
2347: !2244: !
2348: !  Add constraint for this distance to the list.2245: !  Add constraint for this distance to the list.
2349: !2246: !
2350:     DO 2247:     DO 
2351:        READ(LUNIT,*,END=531)  J2, J32248:        READ(LUNIT,*,END=531)  J2, J3, DUMMY
2352: !2249: !
2353: ! Forbid constraints corresponding to atoms distant in sequence. Set INTCONSEP to number of sites to 2250: ! Forbid constraints corresponding to atoms distant in sequence. Set INTCONSEP to number of sites to 
2354: ! turn this off2251: ! turn this off
2355: !2252: !
2356:        IF (J3-J2.GT.INTCONSEP) CYCLE 2253:        IF (J3-J2.GT.INTCONSEP) CYCLE 
2357:        IF (INTFROZEN(J2).AND.INTFROZEN(J3)) CYCLE ! no constraints between intfrozen atoms2254:        IF (INTFROZEN(J2).AND.INTFROZEN(J3)) CYCLE ! no constraints between intfrozen atoms
2358:        NCONSTRAINT=NCONSTRAINT+12255:        NCONSTRAINT=NCONSTRAINT+1
2359:        WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for atoms ',J2,J3,'  total=',NCONSTRAINT2256:        WRITE(MYUNIT,'(A,2I6,A,I6)') 'intlbfgs> Adding constraint for atoms ',J2,J3,'  total=',NCONSTRAINT
2360:        DS=SQRT((LXYZ(3*(J2-1)+1)-LXYZ(3*(J3-1)+1))**2 &2257:        DS=SQRT((LXYZ(3*(J2-1)+1)-LXYZ(3*(J3-1)+1))**2 &
2361:   &           +(LXYZ(3*(J2-1)+2)-LXYZ(3*(J3-1)+2))**2 &2258:   &           +(LXYZ(3*(J2-1)+2)-LXYZ(3*(J3-1)+2))**2 &
2363:        DF=SQRT((LXYZ(3*NATOMS+3*(J2-1)+1)-LXYZ(3*NATOMS+3*(J3-1)+1))**2 &2260:        DF=SQRT((LXYZ(3*NATOMS+3*(J2-1)+1)-LXYZ(3*NATOMS+3*(J3-1)+1))**2 &
2364:   &           +(LXYZ(3*NATOMS+3*(J2-1)+2)-LXYZ(3*NATOMS+3*(J3-1)+2))**2 &2261:   &           +(LXYZ(3*NATOMS+3*(J2-1)+2)-LXYZ(3*NATOMS+3*(J3-1)+2))**2 &
2365:   &           +(LXYZ(3*NATOMS+3*(J2-1)+3)-LXYZ(3*NATOMS+3*(J3-1)+3))**2) 2262:   &           +(LXYZ(3*NATOMS+3*(J2-1)+3)-LXYZ(3*NATOMS+3*(J3-1)+3))**2) 
2366:        IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE2263:        IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE
2367:        CONI(NCONSTRAINT)=J22264:        CONI(NCONSTRAINT)=J2
2368:        CONJ(NCONSTRAINT)=J32265:        CONJ(NCONSTRAINT)=J3
2369:        CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D02266:        CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D0
2370:        CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D02267:        CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D0
2371:        IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)2268:        IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)
2372:        IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)2269:        IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)
2373:        WRITE(MYUNIT,'(A,2I6,A,2F12.2,A,F12.4,A,I8)') ' intlbfgs> constrain distance for atoms ',CONI(NCONSTRAINT), &2270:        WRITE(MYUNIT,'(A,2I6,A,3F12.2,A,F12.4,A,I8)') ' intlbfgs> constrain distance for atoms ',CONI(NCONSTRAINT), &
2374:   &                 CONJ(NCONSTRAINT),' values are ',DS,DF,' fraction=',2*ABS(DS-DF)/(DS+DF), &2271:   &                 CONJ(NCONSTRAINT),' values are ',DS,DF,DUMMY,' fraction=',2*ABS(DS-DF)/(DS+DF), &
2375:   &                ' # constraints=',NCONSTRAINT2272:   &                ' # constraints=',NCONSTRAINT
2376:     ENDDO2273:     ENDDO
2377: 531 CONTINUE2274: 531 CONTINUE
2378:     WRITE(MYUNIT,'(A,I6,2(A,F15.5))') ' intlbfgs> Total distance constraints=',NCONSTRAINT, &2275:     WRITE(MYUNIT,'(A,I6,2(A,F15.5))') ' intlbfgs> Total distance constraints=',NCONSTRAINT, &
2379:   &                               ' shortest=',MINCONDIST,' longest=',MAXCONDIST2276:   &                               ' shortest=',MINCONDIST,' longest=',MAXCONDIST
2380: 2277: 
 2278: 
2381: ELSE IF (NCONGEOM.LT.2) THEN 2279: ELSE IF (NCONGEOM.LT.2) THEN 
2382:    DO J2=1,NATOMS2280:    DO J2=1,NATOMS
2383:       DO J3=J2+1,NATOMS2281:       DO J3=J2+1,NATOMS
2384: 2282: 
2385:          IF (J3-J2.GT.INTCONSEP) CYCLE ! forbid constraints corresponding to atoms distant in sequence2283:          IF (J3-J2.GT.INTCONSEP) CYCLE ! forbid constraints corresponding to atoms distant in sequence
2386:          IF (INTFROZEN(J2).AND.INTFROZEN(J3)) CYCLE ! no constraints between intfrozen atoms2284:          IF (INTFROZEN(J2).AND.INTFROZEN(J3)) CYCLE ! no constraints between intfrozen atoms
2387:          DS=SQRT((LXYZ(3*(J2-1)+1)-LXYZ(3*(J3-1)+1))**2 &2285:          DS=SQRT((LXYZ(3*(J2-1)+1)-LXYZ(3*(J3-1)+1))**2 &
2388:   &             +(LXYZ(3*(J2-1)+2)-LXYZ(3*(J3-1)+2))**2 &2286:   &             +(LXYZ(3*(J2-1)+2)-LXYZ(3*(J3-1)+2))**2 &
2389:   &             +(LXYZ(3*(J2-1)+3)-LXYZ(3*(J3-1)+3))**2) 2287:   &             +(LXYZ(3*(J2-1)+3)-LXYZ(3*(J3-1)+3))**2) 
2390:          IF (DS.GT.INTCONCUT) CYCLE ! don't allow constraints if either endpoint separation is too large DJW2288:          IF (DS.GT.INTCONCUT) CYCLE ! don't allow constraints if either endpoint separation is too large DJW
2450:   &                       ' cutoff=',CONCUT(NCONSTRAINT),' constraints=',NCONSTRAINT2348:   &                       ' cutoff=',CONCUT(NCONSTRAINT),' constraints=',NCONSTRAINT
2451: 753      CONTINUE2349: 753      CONTINUE
2452:       ENDDO2350:       ENDDO
2453:    ENDDO2351:    ENDDO
2454:    CONIFIX(1:NCONSTRAINT)=CONI(1:NCONSTRAINT)2352:    CONIFIX(1:NCONSTRAINT)=CONI(1:NCONSTRAINT)
2455:    CONJFIX(1:NCONSTRAINT)=CONJ(1:NCONSTRAINT)2353:    CONJFIX(1:NCONSTRAINT)=CONJ(1:NCONSTRAINT)
2456:    CONDISTREFFIX(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)2354:    CONDISTREFFIX(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)
2457:    CONCUTFIX(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT)2355:    CONCUTFIX(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT)
2458:    NCONSTRAINTFIX=NCONSTRAINT2356:    NCONSTRAINTFIX=NCONSTRAINT
2459: ENDIF2357: ENDIF
2460:  
2461: write(myunit,'(A,2I6)') 'QCIADDREP,QCIBONDS=',QCIADDREP,QCIBONDS 
2462:  
2463: IF (QCIADDREP.GT.0) THEN 
2464:    DMIN=1.0D100 
2465:    DO J2=1,QCIBONDS 
2466: ! 
2467: ! end point 1 
2468: ! 
2469:       NI1=3*(CONI(J2)-1) 
2470:       NJ1=3*(CONJ(J2)-1) 
2471:       DO J3=J2+1,QCIBONDS 
2472:          IF (CONI(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom 
2473:          IF (CONI(J3).EQ.CONJ(J2)) CYCLE ! no extra terms for bonds with a common atom 
2474:          IF (CONJ(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom 
2475:          IF (CONJ(J3).EQ.CONJ(J2)) CYCLE ! no extra terms for bonds with a common atom 
2476: ! 
2477: ! end point 1 
2478: ! 
2479:          NI2=3*(CONI(J3)-1) 
2480:          NJ2=3*(CONJ(J3)-1) 
2481:          DO J4=1,QCIADDREP 
2482:             X1=(J4*LXYZ(NI1+1)+(QCIADDREP+1-J4)*LXYZ(NJ1+1))/(QCIADDREP+1.0D0) 
2483:             Y1=(J4*LXYZ(NI1+2)+(QCIADDREP+1-J4)*LXYZ(NJ1+2))/(QCIADDREP+1.0D0) 
2484:             Z1=(J4*LXYZ(NI1+3)+(QCIADDREP+1-J4)*LXYZ(NJ1+3))/(QCIADDREP+1.0D0) 
2485:             DO J5=1,QCIADDREP 
2486:                X2=(J5*LXYZ(NI2+1)+(QCIADDREP+1-J5)*LXYZ(NJ2+1))/(QCIADDREP+1.0D0) 
2487:                Y2=(J5*LXYZ(NI2+2)+(QCIADDREP+1-J5)*LXYZ(NJ2+2))/(QCIADDREP+1.0D0) 
2488:                Z2=(J5*LXYZ(NI2+3)+(QCIADDREP+1-J5)*LXYZ(NJ2+3))/(QCIADDREP+1.0D0) 
2489:                D2=SQRT((X1-X2)**2+(Y1-Y2)**2+(Z1-Z2)**2) 
2490:                IF (D2.LT.DMIN) DMIN=D2 
2491: !              WRITE(MYUNIT,'(A,2I6,A,4I6,A,2I6,A,F20.10)') 'intlbfgs> start constraints ',J2,J3,' atoms ', & 
2492: ! &                                CONI(J2),CONJ(J2),CONI(J3),CONJ(J3),' J4,J5 ',J4,J5,' distance=',D2 
2493:            ENDDO 
2494:          ENDDO 
2495:       ENDDO 
2496: ! 
2497: ! end point 2 
2498: ! 
2499:       NI1=3*(CONI(J2)-1)+3*NATOMS 
2500:       NJ1=3*(CONJ(J2)-1)+3*NATOMS 
2501:       DO J3=J2+1,QCIBONDS 
2502:          IF (CONI(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom 
2503:          IF (CONI(J3).EQ.CONJ(J2)) CYCLE ! no extra terms for bonds with a common atom 
2504:          IF (CONJ(J3).EQ.CONI(J2)) CYCLE ! no extra terms for bonds with a common atom 
2505:          IF (CONJ(J3).EQ.CONJ(J2)) CYCLE ! no extra terms for bonds with a common atom 
2506: ! 
2507: ! end point 2 
2508: ! 
2509:      NI2=3*(CONI(J3)-1) 
2510:          NI2=3*(CONI(J3)-1)+3*NATOMS 
2511:          NJ2=3*(CONJ(J3)-1)+3*NATOMS 
2512:          DO J4=1,QCIADDREP 
2513:             X1=(J4*LXYZ(NI1+1)+(QCIADDREP+1-J4)*LXYZ(NJ1+1))/(QCIADDREP+1.0D0) 
2514:             Y1=(J4*LXYZ(NI1+2)+(QCIADDREP+1-J4)*LXYZ(NJ1+2))/(QCIADDREP+1.0D0) 
2515:             Z1=(J4*LXYZ(NI1+3)+(QCIADDREP+1-J4)*LXYZ(NJ1+3))/(QCIADDREP+1.0D0) 
2516:             DO J5=1,QCIADDREP 
2517:                X2=(J5*LXYZ(NI2+1)+(QCIADDREP+1-J5)*LXYZ(NJ2+1))/(QCIADDREP+1.0D0) 
2518:                Y2=(J5*LXYZ(NI2+2)+(QCIADDREP+1-J5)*LXYZ(NJ2+2))/(QCIADDREP+1.0D0) 
2519:                Z2=(J5*LXYZ(NI2+3)+(QCIADDREP+1-J5)*LXYZ(NJ2+3))/(QCIADDREP+1.0D0) 
2520:                D2=SQRT((X1-X2)**2+(Y1-Y2)**2+(Z1-Z2)**2) 
2521:                IF (D2.LT.DMIN) DMIN=D2 
2522: !              WRITE(MYUNIT,'(A,2I6,A,4I6,A,2I6,A,F20.10)') 'intlbfgs> finish constraints ',J2,J3,' atoms ', & 
2523: ! &                                CONI(J2),CONJ(J2),CONI(J3),CONJ(J3),' J4,J5 ',J4,J5,' distance=',D2 
2524:            ENDDO 
2525:          ENDDO 
2526:       ENDDO 
2527:    ENDDO 
2528:    WRITE(MYUNIT,'(A,F20.10,A,F20.10)') 'intlbfgs> minimum decoration distance=',DMIN,' compared with cutoff ',QCIADDREPCUT 
2529:    QCIADDREPCUT=MIN(DMIN-1.0D-3,QCIADDREPCUT) 
2530:    WRITE(MYUNIT,'(A,F20.10)') 'intlbfgs> cutoff after setup is ',QCIADDREPCUT 
2531: ENDIF 
2532: !2358: !
2533: ! Check that we have a percolating constraint network. If not, increase the tolerance and try again!2359: ! Check that we have a percolating constraint network. If not, increase the tolerance and try again!
2534: ! Calculate minimum number of steps of each atom from number 1 or any frozen atom.2360: ! Calculate minimum number of steps of each atom from number 1 or any frozen atom.
2535: !2361: !
2536: NDIST1(1:NATOMS)=10000002362: NDIST1(1:NATOMS)=1000000
2537: IF (NQCIFREEZE.EQ.0) THEN2363: IF (NQCIFREEZE.EQ.0) THEN
2538:    NDIST1(1)=02364:    NDIST1(1)=0
2539: ELSE2365: ELSE
2540:    DO J1=1,NATOMS2366:    DO J1=1,NATOMS
2541:       IF (INTFROZEN(J1)) NDIST1(J1)=02367:       IF (INTFROZEN(J1)) NDIST1(J1)=0


r30603/keywords.f 2016-06-16 09:30:16.025138392 +0100 r30602/keywords.f 2016-06-16 09:30:17.377156539 +0100
165:       CHFREQBB=1165:       CHFREQBB=1
166:       FTRANS=1166:       FTRANS=1
167:       FROT=1167:       FROT=1
168: 168: 
169: !169: !
170: ! QCI parameters170: ! QCI parameters
171: !171: !
172:          CONDATT=.FALSE.172:          CONDATT=.FALSE.
173:          QCIPOTT=.FALSE.173:          QCIPOTT=.FALSE.
174:          QCIPOT2T=.FALSE.174:          QCIPOT2T=.FALSE.
175:          QCIADDREP=0 
176:          QCIADDREPCUT=1.0D0 
177:          QCIADDREPEPS=1.0D0 
178:          QCINOREPINT=.FALSE. 
179:          MAXNACTIVE=0175:          MAXNACTIVE=0
180:          FREEZETOL=1.0D-3176:          FREEZETOL=1.0D-3
181:          QCIPERMCHECK=.FALSE.177:          QCIPERMCHECK=.FALSE.
182:          QCIPERMCHECKINT=100178:          QCIPERMCHECKINT=100
183:          INTCONSTRAINTT=.FALSE.179:          INTCONSTRAINTT=.FALSE.
184:          INTCONSTRAINTTOL=0.1D0180:          INTCONSTRAINTTOL=0.1D0
185:          INTCONSTRAINTDEL=10.0D0181:          INTCONSTRAINTDEL=10.0D0
186:          INTCONSTRAINTREP=100.0D0182:          INTCONSTRAINTREP=100.0D0
187:          INTCONSTRAINREPCUT=1.7D0183:          INTCONSTRAINREPCUT=1.7D0
188:          INTFREEZET=.FALSE.184:          INTFREEZET=.FALSE.
1962:          IF (NITEMS.EQ.2) THEN1958:          IF (NITEMS.EQ.2) THEN
1963:             CALL READF(CENTX) 1959:             CALL READF(CENTX) 
1964:          ELSE IF (NITEMS.EQ.3) THEN1960:          ELSE IF (NITEMS.EQ.3) THEN
1965:             CALL READF(CENTX)1961:             CALL READF(CENTX)
1966:             CALL READF(CENTY)1962:             CALL READF(CENTY)
1967:          ELSE IF (NITEMS.EQ.4) THEN1963:          ELSE IF (NITEMS.EQ.4) THEN
1968:             CALL READF(CENTX)1964:             CALL READF(CENTX)
1969:             CALL READF(CENTY)1965:             CALL READF(CENTY)
1970:             CALL READF(CENTZ)1966:             CALL READF(CENTZ)
1971:          ENDIF1967:          ENDIF
1972:       ELSE IF (WORD.EQ.'QCIADDREP') THEN 
1973:          CALL READI(QCIADDREP) 
1974:          CALL READF(QCIADDREPEPS) 
1975:          CALL READF(QCIADDREPCUT) 
1976:          WRITE(MYUNIT,'(A,I6,A)') 'keywords> Adding ',QCIADDREP,' repulsive sites along constraints' 
1977:       ELSE IF (WORD.EQ.'QCINOREPINT') THEN 
1978:          QCINOREPINT=.TRUE. 
1979:       ELSE IF (WORD.EQ.'QCIPOT') THEN1968:       ELSE IF (WORD.EQ.'QCIPOT') THEN
1980:          QCIPOTT=.TRUE.1969:          QCIPOTT=.TRUE.
1981:       ELSE IF (WORD.EQ.'QCIPOT2') THEN1970:       ELSE IF (WORD.EQ.'QCIPOT2') THEN
1982:          QCIPOTT=.TRUE.1971:          QCIPOTT=.TRUE.
1983:          QCIPOT2T=.TRUE.1972:          QCIPOT2T=.TRUE.
1984: 1973: 
1985: !       csw34> QUCENTRE moves the centre of coordinates to (0,0,0)1974: !       csw34> QUCENTRE moves the centre of coordinates to (0,0,0)
1986: !       before each step is taken. 1975: !       before each step is taken. 
1987:       ELSE IF (WORD.EQ.'QUCENTRE') THEN1976:       ELSE IF (WORD.EQ.'QUCENTRE') THEN
1988:          QUCENTRE=.TRUE.1977:          QUCENTRE=.TRUE.


r30603/make_conpot.f90 2016-06-16 09:30:16.225141072 +0100 r30602/make_conpot.f90 2016-06-16 09:30:17.577159223 +0100
196: ! The rest of this code is for initial setup. It isn't needed if CONDATT is true.196: ! The rest of this code is for initial setup. It isn't needed if CONDATT is true.
197: !197: !
198: REPCON=-INTCONSTRAINTREP/INTCONSTRAINREPCUT198: REPCON=-INTCONSTRAINTREP/INTCONSTRAINREPCUT
199: IF (ALLOCATED(CONDISTREFLOCAL)) DEALLOCATE(CONDISTREFLOCAL)199: IF (ALLOCATED(CONDISTREFLOCAL)) DEALLOCATE(CONDISTREFLOCAL)
200: ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT))200: ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT))
201: CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)201: CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)
202: IF (ALLOCATED(CONCUTLOCAL)) DEALLOCATE(CONCUTLOCAL)202: IF (ALLOCATED(CONCUTLOCAL)) DEALLOCATE(CONCUTLOCAL)
203: ALLOCATE(CONCUTLOCAL(NCONSTRAINT))203: ALLOCATE(CONCUTLOCAL(NCONSTRAINT))
204: CONCUTLOCAL(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT)204: CONCUTLOCAL(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT)
205: NREPULSIVEFIX=0205: NREPULSIVEFIX=0
206: IF (INTCONSTRAINTREP.EQ.0.0D0) GOTO 963 
207: !206: !
208: ! Add repulsions to non-constrained atoms.207: ! Add repulsions to non-constrained atoms.
209: ! Note that we do not limit the number of constraints per site in this208: ! Note that we do not limit the number of constraints per site in this
210: ! routine, unlike NEB/lbfgs.f90, where the result will depend on the209: ! routine, unlike NEB/lbfgs.f90, where the result will depend on the
211: ! order in which the constraints are turned on. 210: ! order in which the constraints are turned on. 
212: !211: !
213: NDUMMY=1212: NDUMMY=1
214: DO J1=1,NATOMS213: DO J1=1,NATOMS
215: !214: !
216: ! Make a list of repelling atoms here and then use it215: ! Make a list of repelling atoms here and then use it
261:       REPI(NREPULSIVEFIX)=J1260:       REPI(NREPULSIVEFIX)=J1
262:       REPJ(NREPULSIVEFIX)=J2261:       REPJ(NREPULSIVEFIX)=J2
263: !262: !
264: ! Use the minimum of the end point distances and INTCONSTRAINREPCUT for each contact.263: ! Use the minimum of the end point distances and INTCONSTRAINREPCUT for each contact.
265: !264: !
266:       REPCUT(NREPULSIVEFIX)=MIN(DMIN-1.0D-3,INTCONSTRAINREPCUT)265:       REPCUT(NREPULSIVEFIX)=MIN(DMIN-1.0D-3,INTCONSTRAINREPCUT)
267:       IF (DEBUG) WRITE(MYUNIT, '(A,I6,A,I6,A,F15.5,A,I10)') ' make_conpot> Adding repulsion for atom ',J1, &266:       IF (DEBUG) WRITE(MYUNIT, '(A,I6,A,I6,A,F15.5,A,I10)') ' make_conpot> Adding repulsion for atom ',J1, &
268:   &              ' with atom ',J2,' cutoff=',DMIN,' # repulsions ',NREPULSIVEFIX267:   &              ' with atom ',J2,' cutoff=',DMIN,' # repulsions ',NREPULSIVEFIX
269:    ENDDO rep2268:    ENDDO rep2
270: ENDDO269: ENDDO
271: 963 CONTINUE 
272: NREPULSIVE=NREPULSIVEFIX270: NREPULSIVE=NREPULSIVEFIX
273: 271: 
274: WRITE(MYUNIT, '(A,2I10,A,G20.10)') ' make_conpot> Total number of constraints and repulsions=', &272: WRITE(MYUNIT, '(A,2I10,A,G20.10)') ' make_conpot> Total number of constraints and repulsions=', &
275:   &   NCONSTRAINT,NREPULSIVE,' for tolerance parameter ',LINTCONSTRAINTTOL273:   &   NCONSTRAINT,NREPULSIVE,' for tolerance parameter ',LINTCONSTRAINTTOL
276: 274: 
277: IF (ALLOCATED(CONACTIVE)) DEALLOCATE(CONACTIVE)275: IF (ALLOCATED(CONACTIVE)) DEALLOCATE(CONACTIVE)
278: ALLOCATE(CONACTIVE(NCONSTRAINT))276: ALLOCATE(CONACTIVE(NCONSTRAINT))
279: CONACTIVE(1:NCONSTRAINT)=.TRUE. 277: CONACTIVE(1:NCONSTRAINT)=.TRUE. 
280: !278: !
281: ! congrad routines actually use NREPI, NREPJ, etc., so we must assign these.279: ! congrad routines actually use NREPI, NREPJ, etc., so we must assign these.


r30603/qcipot.f90 2016-06-16 09:30:16.445144036 +0100 r30602/qcipot.f90 2016-06-16 09:30:17.781161994 +0100
 99:          GGG(3*(J2-1)+1)=0.0D0 99:          GGG(3*(J2-1)+1)=0.0D0
100:          GGG(3*(J2-1)+2)=0.0D0100:          GGG(3*(J2-1)+2)=0.0D0
101:          GGG(3*(J2-1)+3)=0.0D0101:          GGG(3*(J2-1)+3)=0.0D0
102:       ENDIF102:       ENDIF
103:    ENDDO103:    ENDDO
104: ENDIF104: ENDIF
105: 105: 
106: END SUBROUTINE QCIPOT106: END SUBROUTINE QCIPOT
107: 107: 
108: SUBROUTINE CHECKREPQCIPOT(XYZ,NOPT,NNSTART,NSTART)108: SUBROUTINE CHECKREPQCIPOT(XYZ,NOPT,NNSTART,NSTART)
109: USE COMMONS,ONLY : NREPI, NREPJ, NREPCUT, NNREPULSIVE, NREPULSIVE, REPI, REPJ, REPCUT, CHECKREPCUTOFF, DEBUG, MYUNIT, intconstraintrep109: USE COMMONS,ONLY : NREPI, NREPJ, NREPCUT, NNREPULSIVE, NREPULSIVE, REPI, REPJ, REPCUT, CHECKREPCUTOFF, DEBUG, MYUNIT
110: USE PORFUNCS110: USE PORFUNCS
111: IMPLICIT NONE111: IMPLICIT NONE
112: INTEGER JJ, NOPT, NI, NJ, NNSTART, NSTART112: INTEGER JJ, NOPT, NI, NJ, NNSTART, NSTART
113: DOUBLE PRECISION LDIST, XYZ(NOPT),COMPARE113: DOUBLE PRECISION LDIST, XYZ(NOPT),COMPARE
114: 114: 
115: IF (INTCONSTRAINTREP.EQ.0) THEN 
116:    NNREPULSIVE=0 
117:    RETURN 
118: ENDIF 
119:  
120: NNREPULSIVE=NNSTART115: NNREPULSIVE=NNSTART
121: DO JJ=NSTART,NREPULSIVE116: DO JJ=NSTART,NREPULSIVE
122:    COMPARE=(CHECKREPCUTOFF*REPCUT(JJ))**2117:    COMPARE=(CHECKREPCUTOFF*REPCUT(JJ))**2
123:    NI=REPI(JJ)118:    NI=REPI(JJ)
124:    NJ=REPJ(JJ)119:    NJ=REPJ(JJ)
125:    LDIST=(XYZ(3*(NI-1)+1)-XYZ(3*(NJ-1)+1))**2 &120:    LDIST=(XYZ(3*(NI-1)+1)-XYZ(3*(NJ-1)+1))**2 &
126:   &     +(XYZ(3*(NI-1)+2)-XYZ(3*(NJ-1)+2))**2 &121:   &     +(XYZ(3*(NI-1)+2)-XYZ(3*(NJ-1)+2))**2 &
127:   &     +(XYZ(3*(NI-1)+3)-XYZ(3*(NJ-1)+3))**2122:   &     +(XYZ(3*(NI-1)+3)-XYZ(3*(NJ-1)+3))**2
128:    IF (LDIST.LT.COMPARE) THEN123:    IF (LDIST.LT.COMPARE) THEN
129:       NNREPULSIVE=NNREPULSIVE+1124:       NNREPULSIVE=NNREPULSIVE+1


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0