hdiff output

r31263/atom_label_swaps.f90 2016-10-12 11:30:23.820781990 +0100 r31262/atom_label_swaps.f90 2016-10-12 11:30:29.676860577 +0100
146:   INVATOMLISTS(:,:) = INVATOMLISTS_MIN(:,:)146:   INVATOMLISTS(:,:) = INVATOMLISTS_MIN(:,:)
147:   COORDS(1:3*NATOMS, NP) = XMIN(1:3*NATOMS)147:   COORDS(1:3*NATOMS, NP) = XMIN(1:3*NATOMS)
148:   POTEL = EMIN148:   POTEL = EMIN
149:   !149:   !
150:   WRITE(MYUNIT,'(A,G20.10,A,F8.5)') &150:   WRITE(MYUNIT,'(A,G20.10,A,F8.5)') &
151:        'uswaps> Final E= ', EMIN,' T= ', TEMP151:        'uswaps> Final E= ', EMIN,' T= ', TEMP
152:   !152:   !
153:   RETURN153:   RETURN
154:   !154:   !
155: END SUBROUTINE USWAPS155: END SUBROUTINE USWAPS
156:  
157: SUBROUTINE SWAP_COORDS_V2(X,I,J) 
158:   ! 
159:   USE COMMONS, ONLY : NATOMS 
160:   ! 
161:   IMPLICIT NONE 
162:   ! 
163:   INTEGER, INTENT(IN) :: I,J 
164:   DOUBLE PRECISION, INTENT(INOUT) :: X(1:3*NATOMS) 
165:   INTEGER :: K,L,M 
166:   DOUBLE PRECISION :: DUMMY 
167:   ! 
168:   K = 3*(I-1) 
169:   L = 3*(J-1) 
170:   DO M = 1,3 
171:      DUMMY = X(K+M) 
172:      X(K+M) = X(L+M) 
173:      X(L+M) = DUMMY 
174:   ENDDO 
175:   ! 
176: END SUBROUTINE SWAP_COORDS_V2 


r31263/commons.f90 2016-10-12 11:30:24.348789174 +0100 r31262/commons.f90 2016-10-12 11:30:30.212867782 +0100
 12: !  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: !  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: !  GNU General Public License for more details. 13: !  GNU General Public License for more details.
 14: ! 14: !
 15: !  You should have received a copy of the GNU General Public License 15: !  You should have received a copy of the GNU General Public License
 16: !  along with this program; if not, write to the Free Software 16: !  along with this program; if not, write to the Free Software
 17: !  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19:       MODULE COMMONS 19:       MODULE COMMONS
 20:       IMPLICIT NONE 20:       IMPLICIT NONE
 21:       SAVE 21:       SAVE
 22:       INTEGER :: NATOMS, NATOMSALLOC, NACCEPT, MAXIT, NTARGETS, NN, MM, NPAR, NPAR_GBH, NPATCH, TIPID, PAHID, NSYMREM=0, NFREEZE, NSAVEINTE, & 22:       INTEGER :: NATOMS, NATOMSALLOC, NACCEPT, MAXIT, NTARGETS, NN, MM, NPAR, NPATCH, TIPID, PAHID, NSYMREM=0, NFREEZE, NSAVEINTE, &
 23:      &        NINTV, NSEED, NVEC, NS, NSSTOP, MAXIT2, NSAVE, SHELLMOVEMAX, HPTGRP, NFIX, CURRENTIMP,& 23:      &        NINTV, NSEED, NVEC, NS, NSSTOP, MAXIT2, NSAVE, SHELLMOVEMAX, HPTGRP, NFIX, CURRENTIMP,&
 24:      &        MCSTEPS(3), NRUNS, NVAR, MYPOWER, XMUL,  NTAB, NCARBON, NCOOP, NSYMINTERVAL, QDLIMIT, COULN, & 24:      &        MCSTEPS(3), NRUNS, NVAR, MYPOWER, XMUL,  NTAB, NCARBON, NCOOP, NSYMINTERVAL, QDLIMIT, COULN, &
 25:      &        NSUPER, NSACCEPT, NSUPERSTEP, GMODES, QUENCHFRQ, CHNRES, CHECKDID, & 25:      &        NSUPER, NSACCEPT, NSUPERSTEP, GMODES, QUENCHFRQ, CHNRES, CHECKDID, &
 26:      &        NEVL, NEVS, NVECTORS, MUPDATE, LSYS, NHSMOVE, NRENORM, NRELAX, NHSRESTART, NRENSTUCK,& 26:      &        NEVL, NEVS, NVECTORS, MUPDATE, LSYS, NHSMOVE, NRENORM, NRELAX, NHSRESTART, NRENSTUCK,&
 27:      &        HBINS, NRBSITES, NRBSITES1, NPS, GATOM, NMSBSAVE, MAXSAVE, NSYMQMAX, EQUIL, TPAHA, ISTEP, & 27:      &        HBINS, NRBSITES, NRBSITES1, NPS, GATOM, NMSBSAVE, MAXSAVE, NSYMQMAX, EQUIL, TPAHA, ISTEP, &
 28:      &        NINTS, NumberOfTrajectories, & 28:      &        NINTS, NumberOfTrajectories, &
 29:      &        MCCycles, MCCyclesEquilibration, TargetWL, NTempPoints, SaveNth, DumpEveryNthQuench, NSpline, Hwindows, &  29:      &        MCCycles, MCCyclesEquilibration, TargetWL, NTempPoints, SaveNth, DumpEveryNthQuench, NSpline, Hwindows, & 
 30:      &        lhbins, sampledbins, MYNODE, MYNQ, NENRPER, DUMPINT, MYUNIT, PRTFRQ, BSPTDUMPFRQ, NPERMGROUP, NMOVABLEATOMS, & 30:      &        lhbins, sampledbins, MYNODE, NENRPER, DUMPINT, MYUNIT, PRTFRQ, BSPTDUMPFRQ, NPERMGROUP, NMOVABLEATOMS, &
 31:      &        CPS, CPF, ACKLANDID, PATOM1, PATOM2, CSMGPINDEX, CSMGUIDEGPINDEX, CSMSTEPS, CSMQUENCHES, CSMMAXIT, & 31:      &        CPS, CPF, ACKLANDID, PATOM1, PATOM2, CSMGPINDEX, CSMGUIDEGPINDEX, CSMSTEPS, CSMQUENCHES, CSMMAXIT, &
 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, &
 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 86:      &                 QCIADDREPCUT, QCIADDREPEPS
 87:  87: 
 88:       LOGICAL DEBUG, TARGET, MORSET, CUTT, SEEDT, CENT, TSALLIST, FREEZECORE, NEWJUMP, RENORM, CAPSID, FREEZE, & 88:       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, & 89:      &        OTPT, LJMFT, STRANDT, PAHT, SWT, MSTRANST, STOCKT, STICKYT, BLNT, MYSDT, FREEZERES, CENTXY, &
 90:      &        MSORIGT, SQUEEZET, PERIODIC, SCT, MSCT, MGUPTAT, RESIZET, TIP, RIGID, CALCQT, MPIT, GBHT, JMT, LJCOULT, SETCENT, & 90:      &        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, & 91:      &        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, & 92:      &        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, & 93:      &        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, & 94:      &        EAMALT, ALGLUET, MGGLUET, GUPTAT, LJATT, FST, DECAY, COOP, FIXBIN, GAUSST, QUENCHDOS, FIXDIHEFLAG, &
 95:      &        FRAUSIT, ANGST, SELFT, STEPOUT, WENZEL, THRESHOLDT, THOMSONT, MULLERBROWNT, CHARMMENERGIES, & 95:      &        FRAUSIT, ANGST, SELFT, STEPOUT, WENZEL, THRESHOLDT, THOMSONT, MULLERBROWNT, CHARMMENERGIES, &
 96:      &        PROJ, RGCL2, TOSI, WELCH, AXTELL, AMBER, FIXIMAGE, BINARY, SHIFTCUT, ARNO, TUNNELT, TWOD, &  96:      &        PROJ, RGCL2, TOSI, WELCH, AXTELL, AMBER, FIXIMAGE, BINARY, SHIFTCUT, ARNO, TUNNELT, TWOD, & 
 97:      &        BLJCLUSTER, BLJCLUSTER_NOCUT, COMPRESST, FIX, FIXT, BFGS, LBFGST, DBRENTT, DZTEST, FNI, FAL, CPMD, TNT, ZETT1, & 97:      &        BLJCLUSTER, BLJCLUSTER_NOCUT, COMPRESST, FIX, FIXT, BFGS, LBFGST, DBRENTT, DZTEST, FNI, FAL, CPMD, TNT, ZETT1, &
 98:      &        ZETT2, RESTART, CONJG, NEWRESTART, AVOID, NATBT, DIFFRACTT, CHRMMT, INTMINT, LB2T, &  98:      &        ZETT2, RESTART, CONJG, NEWRESTART, AVOID, NATBT, DIFFRACTT, CHRMMT, INTMINT, LB2T, & 
 99:      &        PTMC, BINSTRUCTURES, PROGRESS, MODEL1T, NEWRESTART_MD, CHANGE_TEMP, NOCISTRANS, CHECKCHIRALITY, & 99:      &        PTMC, BINSTRUCTURES, PROGRESS, MODEL1T, NEWRESTART_MD, CHANGE_TEMP, NOCISTRANS, CHECKCHIRALITY, &
100:      &        GBT, GBDT, GBDPT, GEMT, LINRODT, RADIFT, CAPBINT, DBPT, DBPTDT, DMBLMT, DMBLPYT, EFIELDT, PAHAT, STOCKAAT, MORSEDPT, &100:      &        GBT, GBDT, GBDPT, GEMT, LINRODT, RADIFT, CAPBINT, DBPT, DBPTDT, DMBLMT, DMBLPYT, EFIELDT, PAHAT, STOCKAAT, MORSEDPT, &
295:       INTEGER PTEX_INDEP_NSTEPS295:       INTEGER PTEX_INDEP_NSTEPS
296:       INTEGER PTEX_INDEP_UNIT296:       INTEGER PTEX_INDEP_UNIT
297: !js850> 297: !js850> 
298:       INTEGER NPCALL ! the number of potential calls298:       INTEGER NPCALL ! the number of potential calls
299: !js850> for keyword max_npcall299: !js850> for keyword max_npcall
300:       LOGICAL MAX_NPCALLT300:       LOGICAL MAX_NPCALLT
301:       INTEGER MAX_NPCALL301:       INTEGER MAX_NPCALL
302: 302: 
303: !ds656> Binary Gupta keywords, parameters, variables...303: !ds656> Binary Gupta keywords, parameters, variables...
304:       LOGICAL :: BGUPTAT, BGUPTATAB, BGUPTATBB304:       LOGICAL :: BGUPTAT, BGUPTATAB, BGUPTATBB
 305: !ds656> Binary swap keywords and parameters
 306:       LOGICAL :: BASWAP, BASWAPTEST
 307:       INTEGER :: BASWAP_NMAX, BASWAP_NWAIT
 308:       DOUBLE PRECISION :: BASWAP_FRAC !, BASWAP_TEMP no longer used 
305: !ds656> Sequence of label swaps (a.k.a. block of basin-hopping with309: !ds656> Sequence of label swaps (a.k.a. block of basin-hopping with
306: !       exchange moves only)310: !       exchange moves only)
307:       LOGICAL :: LSWAPST, FIXLABELS311:       LOGICAL :: LSWAPST, FIXLABELS
308:       INTEGER :: LSWAPS_N1, LSWAPS_N2, LSWAPS_NUP312:       INTEGER :: LSWAPS_N1, LSWAPS_N2, LSWAPS_NUP
309:       DOUBLE PRECISION :: LSWAPS_TEMP, LSWAPS_TFACTOR313:       DOUBLE PRECISION :: LSWAPS_TEMP, LSWAPS_TFACTOR
310: !ds656> Sequence of label flips (a.k.a. block of semi-grand-314: !ds656> Sequence of label flips (a.k.a. block of semi-grand-
311: !       canonical basin-hopping with flip moves only)315: !       canonical basin-hopping with flip moves only)
312:       LOGICAL :: LFLIPST, LFLIPS_RESET316:       LOGICAL :: LFLIPST, LFLIPS_RESET
313:       INTEGER :: LFLIPS_N1, LFLIPS_N2, LFLIPS_NUP317:       INTEGER :: LFLIPS_N1, LFLIPS_N2, LFLIPS_NUP
314:       DOUBLE PRECISION :: LFLIPS_TEMP, LFLIPS_TFACTOR318:       DOUBLE PRECISION :: LFLIPS_TEMP, LFLIPS_TFACTOR


r31263/countatoms.f90 2016-10-12 11:30:24.604792496 +0100 r31262/countatoms.f90 2016-10-12 11:30:30.476871314 +0100
 21:  21: 
 22: MODULE NOA 22: MODULE NOA
 23:   USE MODAMBER9 , ONLY : INPCRD 23:   USE MODAMBER9 , ONLY : INPCRD
 24:   USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_SETUP 24:   USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_SETUP
 25:   IMPLICIT NONE 25:   IMPLICIT NONE
 26:   SAVE 26:   SAVE
 27:   INTEGER :: NUMBER_OF_ATOMS 27:   INTEGER :: NUMBER_OF_ATOMS
 28:    28:   
 29: CONTAINS 29: CONTAINS
 30:    30:   
 31:   SUBROUTINE COUNTATOMS(MYUNIT, NPAR, GCBHT, GCMU, GCNATOMS, GBHT) 31:   SUBROUTINE COUNTATOMS(MYUNIT, NPAR, GCBHT, GCMU, GCNATOMS)
 32:      32:     
 33:     IMPLICIT NONE 33:     IMPLICIT NONE
 34:     INTEGER, INTENT(IN) :: NPAR 34:     INTEGER, INTENT(IN) :: NPAR
 35:     INTEGER :: EOF,NRES,SEQ(500),I_RES,NOGLY,GLY, MYUNIT,J 35:     INTEGER :: EOF,NRES,SEQ(500),I_RES,NOGLY,GLY, MYUNIT,J
 36:     LOGICAL :: YESNO, YESNOA, YESNOC, YESNOAMH, YESNOA9 36:     LOGICAL :: YESNO, YESNOA, YESNOC, YESNOAMH, YESNOA9
 37:     CHARACTER(LEN=5) TARFL 37:     CHARACTER(LEN=5) TARFL
 38:     CHARACTER(LEN=10)  CHECK 38:     CHARACTER(LEN=10)  CHECK
 39:     CHARACTER(LEN=80) MYLINE,INPCRD1 39:     CHARACTER(LEN=80) MYLINE,INPCRD1
 40:     LOGICAL GCBHT, END, GBHT 40:     LOGICAL GCBHT, END
 41:     DOUBLE PRECISION GCMU 41:     DOUBLE PRECISION GCMU
 42:     INTEGER GCNATOMS, LUNIT, GETUNIT 42:     INTEGER GCNATOMS, LUNIT, GETUNIT
 43:     CHARACTER(LEN=16) WORD 43:     CHARACTER(LEN=16) WORD
 44:      44:     
 45: ! ds656> First scan 'data' for GCBH and GBH, since the latter  
 46: !        affects atom counting. 
 47:  
 48:     LUNIT=GETUNIT() 
 49:     OPEN(UNIT=LUNIT,FILE='data',STATUS='OLD') 
 50:      
 51: 190 CALL INPUT(END,LUNIT) 
 52:     IF (.NOT. END) THEN 
 53:        CALL READU(WORD) 
 54:     ENDIF 
 55:     IF (END .OR. WORD .EQ. 'STOP') GOTO 191 
 56: ! 
 57: ! grand canonical basin-hopping or GBH: 
 58: ! 
 59:     IF (WORD.EQ.'GCBH') THEN 
 60:        GCBHT=.TRUE. 
 61:        CALL READF(GCMU) ! chemical potential 
 62:        CALL READI(GCNATOMS) ! maximum number of atoms 
 63:        WRITE(MYUNIT,'(A,G20.10,A,I6)') & 
 64:             'countatoms> Grand canonical basin-hopping, chemical potential=', & 
 65:             GCMU,' maximum atoms=',GCNATOMS 
 66:        !GOTO 191 ! removed by ds656 to allow complete scan of file 
 67:     ELSE IF (WORD.EQ.'GBH') THEN 
 68:        GBHT=.TRUE. 
 69:        WRITE(MYUNIT,'(A)') 'countatoms> Detected keyword GBH.' 
 70:     ENDIF 
 71:      
 72:     GOTO 190 
 73: 191 CONTINUE 
 74:     CLOSE(LUNIT) 
 75:     ! 
 76:     ! ds656> Now start the atom counting... 
 77:     ! 
 78:     YESNO=.FALSE. 45:     YESNO=.FALSE.
 79:     YESNOA=.FALSE. 46:     YESNOA=.FALSE.
 80:     YESNOAMH=.FALSE. 47:     YESNOAMH=.FALSE.
 81:     YESNOA9=.FALSE. 48:     YESNOA9=.FALSE.
 82:     ! 49:     !
 83:     ! If the current working directory contains more than one of  50:     ! If the current working directory contains more than one of 
 84:     ! these files then the precedence is coords, then input.crd,  51:     ! these files then the precedence is coords, then input.crd, 
 85:     ! then coords.amber  52:     ! then coords.amber 
 86:     ! OPTIM does this a bit better by calling getparams first to  53:     ! OPTIM does this a bit better by calling getparams first to 
 87:     ! see if we are actually doing AMBER or CHARMM.  54:     ! see if we are actually doing AMBER or CHARMM. 
 99:        DO 66:        DO
100:           READ(7,*,IOSTAT=EOF) 67:           READ(7,*,IOSTAT=EOF)
101:           IF (EOF==0) THEN 68:           IF (EOF==0) THEN
102:              NUMBER_OF_ATOMS = NUMBER_OF_ATOMS + 1 69:              NUMBER_OF_ATOMS = NUMBER_OF_ATOMS + 1
103:           ELSE 70:           ELSE
104:              EXIT 71:              EXIT
105:           ENDIF 72:           ENDIF
106:        ENDDO 73:        ENDDO
107:        !js850> if a parallel run, then file coords is NPAR*NATOMS long. 74:        !js850> if a parallel run, then file coords is NPAR*NATOMS long.
108:        !       NPAR must be known at the time of this call 75:        !       NPAR must be known at the time of this call
109:        IF ( NPAR .GT. 1 .AND. .NOT. GBHT) THEN  76:        IF ( NPAR .GT. 1 ) THEN 
110:           IF (MOD(NUMBER_OF_ATOMS,NPAR).NE.0) THEN 77:           IF (MOD(NUMBER_OF_ATOMS,NPAR).NE.0) THEN
111:              WRITE(MYUNIT,'(A,I8,A,I8)') & 78:              WRITE(MYUNIT,'(A,I8,A,I8)') &
112:                   'Number of atoms in coords file=', & 79:                   'Number of atoms in coords file=', &
113:                   NUMBER_OF_ATOMS, & 80:                   NUMBER_OF_ATOMS, &
114:                   ' is not divisible by number of runs=',NPAR 81:                   ' is not divisible by number of runs=',NPAR
115:              STOP 82:              STOP
116:           ENDIF 83:           ENDIF
117:           NUMBER_OF_ATOMS=NUMBER_OF_ATOMS/NPAR 84:           NUMBER_OF_ATOMS=NUMBER_OF_ATOMS/NPAR
118:        ENDIF 85:        ENDIF
119:     ELSEIF (YESNOAMH) THEN 86:     ELSEIF (YESNOAMH) THEN
197:                &file to determine number of atoms. Also could not be guessed from &164:                &file to determine number of atoms. Also could not be guessed from &
198:                &data file (only works for DMACRYS)'165:                &data file (only works for DMACRYS)'
199:           STOP166:           STOP
200:        endif167:        endif
201:     ENDIF168:     ENDIF
202:     169:     
203:     CLOSE(7)170:     CLOSE(7)
204:     171:     
205:     ! ds656> Determine the number of species from 'data'172:     ! ds656> Determine the number of species from 'data'
206:     CALL COUNT_SPECIES_FROM_DATAFILE()173:     CALL COUNT_SPECIES_FROM_DATAFILE()
 174: 
 175:     LUNIT=GETUNIT()
 176:     OPEN(UNIT=LUNIT,FILE='data',STATUS='OLD')
 177:     
 178: 190 CALL INPUT(END,LUNIT)
 179:     IF (.NOT. END) THEN
 180:        CALL READU(WORD)
 181:     ENDIF
 182:     IF (END .OR. WORD .EQ. 'STOP') GOTO 191
 183: !
 184: ! Grand canonical basin-hopping.
 185: !
 186:     IF (WORD.EQ.'GCBH') THEN
 187:        GCBHT=.TRUE.
 188:        CALL READF(GCMU) ! chemical potential
 189:        CALL READI(GCNATOMS) ! maximum number of atoms
 190:        WRITE(MYUNIT,'(A,G20.10,A,I6)') &
 191:             'countatoms> Grand canonical basin-hopping, chemical potential=', &
 192:             GCMU,' maximum atoms=',GCNATOMS
 193:        GOTO 191
 194:     ENDIF
 195:     
 196:     GOTO 190
 197: 191 CONTINUE
 198:     CLOSE(LUNIT)
207:     199:     
208:     !print *, Number_of_Atoms, ' atoms in the system'    200:     !print *, Number_of_Atoms, ' atoms in the system'    
209: 201: 
210:     RETURN202:     RETURN
211:     !203:     !
212:   END SUBROUTINE COUNTATOMS204:   END SUBROUTINE COUNTATOMS
213:   !205:   !
214:   SUBROUTINE COUNT_SPECIES_FROM_DATAFILE()206:   SUBROUTINE COUNT_SPECIES_FROM_DATAFILE()
215:     !207:     !
216:     USE COMMONS, ONLY : NSPECIES, NSPECIES_INI208:     USE COMMONS, ONLY : NSPECIES, NSPECIES_INI


r31263/finalio.f90 2016-10-12 11:30:24.860795941 +0100 r31262/finalio.f90 2016-10-12 11:30:30.740874907 +0100
116:   116:   
117:   IF (DMACRYST) THEN117:   IF (DMACRYST) THEN
118:      CALL DMACRYS_DUMP_LOWEST118:      CALL DMACRYS_DUMP_LOWEST
119:   ENDIF119:   ENDIF
120:   120:   
121:   IF(USERPOTT) THEN121:   IF(USERPOTT) THEN
122:      CALL USERPOT_DUMP_LOWEST122:      CALL USERPOT_DUMP_LOWEST
123:   ENDIF123:   ENDIF
124:   124:   
125:   125:   
126:   IF (MPIT.OR.GBHT) THEN126:   IF (MPIT) THEN
127:      WRITE (ISTR, '(I10)') MYNODE+1127:      WRITE (ISTR, '(I10)') MYNODE+1
128:      MYUNIT2=GETUNIT()128:      MYUNIT2=GETUNIT()
129:      MYFILENAME2="lowest."//TRIM(ADJUSTL(ISTR))129:      MYFILENAME2="lowest."//TRIM(ADJUSTL(ISTR))
130:      !        WRITE(MYUNIT,'(A,I6,A)') ' in finalio MYUNIT2,MYFILENAME2=',MYUNIT2,TRIM(ADJUSTL(MYFILENAME2))130:      !        WRITE(MYUNIT,'(A,I6,A)') ' in finalio MYUNIT2,MYFILENAME2=',MYUNIT2,TRIM(ADJUSTL(MYFILENAME2))
131:      OPEN(MYUNIT2,FILE=TRIM(ADJUSTL(MYFILENAME2)), STATUS="UNKNOWN", FORM="FORMATTED")131:      OPEN(MYUNIT2,FILE=TRIM(ADJUSTL(MYFILENAME2)), STATUS="UNKNOWN", FORM="FORMATTED")
132:      IF (CHRMMT) THEN132:      IF (CHRMMT) THEN
133:         DO J1=1,NSAVE133:         DO J1=1,NSAVE
134:            WRITE(DBNUM,*) J1134:            WRITE(DBNUM,*) J1
135:            DBNAME(J1)='dbase.'//TRIM(ADJUSTL(ISTR))//'.'//TRIM(ADJUSTL(DBNUM))135:            DBNAME(J1)='dbase.'//TRIM(ADJUSTL(ISTR))//'.'//TRIM(ADJUSTL(DBNUM))
136:         ENDDO136:         ENDDO
318:         ELSE ! Print stuff for HSA in OPTIM's min.data format318:         ELSE ! Print stuff for HSA in OPTIM's min.data format
319:            !319:            !
320:            NUM_ZERO_EVS=6320:            NUM_ZERO_EVS=6
321:            IF (NATOMS.EQ.2) NUM_ZERO_EVS=5321:            IF (NATOMS.EQ.2) NUM_ZERO_EVS=5
322:            IF (MIEFT) NUM_ZERO_EVS=0322:            IF (MIEFT) NUM_ZERO_EVS=0
323:            !323:            !
324:            IF (ALLOCATED(HESS)) DEALLOCATE(HESS)324:            IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
325:            ALLOCATE(HESS(3*NATOMS,3*NATOMS))325:            ALLOCATE(HESS(3*NATOMS,3*NATOMS))
326:            IF(NSPECIES(0)>1) THEN ! ds656> this also sets ATMASS326:            IF(NSPECIES(0)>1) THEN ! ds656> this also sets ATMASS
327:               CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)327:               CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)
328:               CALL SET_LABELS(QMINT(J1,1:NATOMS),1)328:               CALL SET_LABELS(QMINT(J1,1:NATOMS),MYNODE+1)
329:            ENDIF329:            ENDIF
330:            CALL POTENTIAL(QMINP(J1,1:3*NATOMS), EVALS, EDUMMY, &330:            CALL POTENTIAL(QMINP(J1,1:3*NATOMS), EVALS, EDUMMY, &
331:                 .TRUE., .TRUE.) ! EVALS is dummy for gradient here331:                 .TRUE., .TRUE.) ! EVALS is dummy for gradient here
332:            IF(ABS(EDUMMY-QMIN(J1)) > ECONV) THEN332:            IF(ABS(EDUMMY-QMIN(J1)) > ECONV) THEN
333:               WRITE(MYUNIT,'(A)') &333:               WRITE(MYUNIT,'(A)') &
334:                    'finalio> WARNING: Unexpected energy value! Configuration changed on final requench?'334:                    'finalio> WARNING: Unexpected energy value! Configuration changed on final requench?'
335:            ENDIF335:            ENDIF
336:            CALL MASSWT() ! Weight the Hessian by mass. 336:            CALL MASSWT() ! Weight the Hessian by mass. 
337:            !337:            !
338:            WRITE(MYUNIT, '(A)') &338:            WRITE(MYUNIT, '(A)') &
575:         END DO575:         END DO
576:         576:         
577:      ELSE IF (AMBER.OR.MOLECULART) THEN577:      ELSE IF (AMBER.OR.MOLECULART) THEN
578:         DO J2=1,NATOMS578:         DO J2=1,NATOMS
579:            WRITE(MYUNIT2,'(A,3F20.10)') typech(J2)(1:1),(QMINP(J1,3*(J2-1)+J3),J3=1,3)579:            WRITE(MYUNIT2,'(A,3F20.10)') typech(J2)(1:1),(QMINP(J1,3*(J2-1)+J3),J3=1,3)
580:         ENDDO580:         ENDDO
581:      ELSE IF (AMBER12T) THEN581:      ELSE IF (AMBER12T) THEN
582:         ! Create a string for J1582:         ! Create a string for J1
583:         WRITE(J1_STRING,'(I6)') J1583:         WRITE(J1_STRING,'(I6)') J1
584:         IF (DUMPSTRUCTURES) THEN584:         IF (DUMPSTRUCTURES) THEN
585:            IF (MPIT.OR.GBHT) THEN585:            IF (MPIT) THEN
586:               ! If we're writing output for multiple parallel MPI runs we need to use MYNODE to586:               ! If we're writing output for multiple parallel MPI runs we need to use MYNODE to
587:               ! distinguish the outputs.587:               ! distinguish the outputs.
588:               ! Create a string for the node too.588:               ! Create a string for the node too.
589:               WRITE(MYNODE_STRING,'(I6)') MYNODE589:               WRITE(MYNODE_STRING,'(I6)') MYNODE
590:               CALL AMBER12_WRITE_RESTART(QMINP(J1,:), 'coords.'//TRIM(ADJUSTL(J1_STRING))//&590:               CALL AMBER12_WRITE_RESTART(QMINP(J1,:), 'coords.'//TRIM(ADJUSTL(J1_STRING))//&
591:                    &'.'//TRIM(ADJUSTL(MYNODE_STRING))//&591:                    &'.'//TRIM(ADJUSTL(MYNODE_STRING))//&
592:                    &'.rst', &592:                    &'.rst', &
593:                    &  LEN('coords.'//TRIM(ADJUSTL(J1_STRING))//&593:                    &  LEN('coords.'//TRIM(ADJUSTL(J1_STRING))//&
594:                    &'.'//TRIM(ADJUSTL(MYNODE_STRING))//&594:                    &'.'//TRIM(ADJUSTL(MYNODE_STRING))//&
595:                    &'.rst'))595:                    &'.rst'))
602:            ELSE602:            ELSE
603:               ! Otherwise, just use one output without the MYNODE suffix.603:               ! Otherwise, just use one output without the MYNODE suffix.
604:               CALL AMBER12_WRITE_RESTART(QMINP(J1,:), 'coords.'//TRIM(ADJUSTL(J1_STRING))//&604:               CALL AMBER12_WRITE_RESTART(QMINP(J1,:), 'coords.'//TRIM(ADJUSTL(J1_STRING))//&
605:                    &'.rst', &605:                    &'.rst', &
606:                    & LEN('coords.'//TRIM(ADJUSTL(J1_STRING))//'.rst'))606:                    & LEN('coords.'//TRIM(ADJUSTL(J1_STRING))//'.rst'))
607:               CALL AMBER12_WRITE_PDB(QMINP(J1,:), 'coords.'//TRIM(ADJUSTL(J1_STRING))//&607:               CALL AMBER12_WRITE_PDB(QMINP(J1,:), 'coords.'//TRIM(ADJUSTL(J1_STRING))//&
608:                    &'.pdb', &608:                    &'.pdb', &
609:                    & LEN('coords.'//TRIM(ADJUSTL(J1_STRING))//'.pdb'))609:                    & LEN('coords.'//TRIM(ADJUSTL(J1_STRING))//'.pdb'))
610:            END IF610:            END IF
611:         END IF611:         END IF
612:         IF (MPIT.OR.GBHT) THEN612:         IF (MPIT) THEN
613:            WRITE(MYNODE_STRING,'(I6)') MYNODE613:            WRITE(MYNODE_STRING,'(I6)') MYNODE
614:            ! Pass the name, length of the string and whether or not to write a header.614:            ! Pass the name, length of the string and whether or not to write a header.
615:            CALL AMBER12_WRITE_XYZ(QMINP(J1,:), 'lowest.'//TRIM(ADJUSTL(MYNODE_STRING)), &615:            CALL AMBER12_WRITE_XYZ(QMINP(J1,:), 'lowest.'//TRIM(ADJUSTL(MYNODE_STRING)), &
616:                 LEN('lowest.'//TRIM(ADJUSTL(MYNODE_STRING))), &616:                 LEN('lowest.'//TRIM(ADJUSTL(MYNODE_STRING))), &
617:                 LOGICAL(.FALSE.,KIND=1))617:                 LOGICAL(.FALSE.,KIND=1))
618:         ELSE 618:         ELSE 
619:            ! Pass the name, length of the string and whether or not to write a header.619:            ! Pass the name, length of the string and whether or not to write a header.
620:            CALL AMBER12_WRITE_XYZ(QMINP(J1,:), 'lowest', LEN('lowest'), LOGICAL(.FALSE.,KIND=1))620:            CALL AMBER12_WRITE_XYZ(QMINP(J1,:), 'lowest', LEN('lowest'), LOGICAL(.FALSE.,KIND=1))
621:         END IF621:         END IF
622:      ELSE IF (AMBERT) THEN622:      ELSE IF (AMBERT) THEN
623:         ! sf344> write out coordinates623:         ! sf344> write out coordinates
624:         COORDS1(1:3*NATOMS) = QMINP(J1,1:3*NATOMS)624:         COORDS1(1:3*NATOMS) = QMINP(J1,1:3*NATOMS)
625:         IF (DUMPSTRUCTURES) THEN625:         IF (DUMPSTRUCTURES) THEN
626:            WRITE(J1CHAR2,'(I3)') J1626:            WRITE(J1CHAR2,'(I3)') J1
627:            IF (MPIT.OR.GBHT) THEN627:            IF (MPIT) THEN
628:               CALL AMBERFINALIO(j1,MYUNIT2,MYNODE+1,tempstring,1,COORDS1(1:3*NATOMS))628:               CALL AMBERFINALIO(j1,MYUNIT2,MYNODE+1,tempstring,1,COORDS1(1:3*NATOMS))
629:               WRITE (ISTR, '(I10)') MYNODE+1629:               WRITE (ISTR, '(I10)') MYNODE+1
630:               MYUNIT3=GETUNIT()630:               MYUNIT3=GETUNIT()
631:               MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))//'.'//trim(adjustl(istr))631:               MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))//'.'//trim(adjustl(istr))
632:               OPEN(MYUNIT3,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")632:               OPEN(MYUNIT3,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")
633:            ELSE633:            ELSE
634:               CALL AMBERFINALIO(j1,MYUNIT2,MYNODE+1,tempstring,0,COORDS1(1:3*NATOMS))634:               CALL AMBERFINALIO(j1,MYUNIT2,MYNODE+1,tempstring,0,COORDS1(1:3*NATOMS))
635:               MYUNIT3=GETUNIT()635:               MYUNIT3=GETUNIT()
636:               MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))636:               MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))
637:               OPEN(MYUNIT3,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")637:               OPEN(MYUNIT3,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")
744:         !744:         !
745:         ! this writes 'lowest' in xyz (Xmakemol) format745:         ! this writes 'lowest' in xyz (Xmakemol) format
746:         !746:         !
747:         WRITE(MYUNIT,'(A,I6,A)') ' in finalio BLN block MYUNIT2=',MYUNIT2747:         WRITE(MYUNIT,'(A,I6,A)') ' in finalio BLN block MYUNIT2=',MYUNIT2
748:         DO J2=1,NATOMS748:         DO J2=1,NATOMS
749:            WRITE(MYUNIT2,'(2A1,1X,3F20.10)') BEADLETTER(J2),'L',(QMINP(J1,3*(J2-1)+J3),J3=1,3)749:            WRITE(MYUNIT2,'(2A1,1X,3F20.10)') BEADLETTER(J2),'L',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
750:         ENDDO750:         ENDDO
751:      ELSE IF(QALCST.OR.MULTIPERMT.OR.MLJT.OR.MGUPTAT.OR.MSCT) THEN751:      ELSE IF(QALCST.OR.MULTIPERMT.OR.MLJT.OR.MGUPTAT.OR.MSCT) THEN
752:         !752:         !
753:         CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)753:         CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)
754:         CALL SET_LABELS(QMINT(J1,1:NATOMS),1)754:         CALL SET_LABELS(QMINT(J1,1:NATOMS),MYNODE+1)
755:         IF(SPECMASST) THEN 755:         IF(SPECMASST) THEN 
756:            EDUMMY = SUM(ATMASS)756:            EDUMMY = SUM(ATMASS)
757:            WRITE(MYUNIT,'(A,I6,A,E20.10)') &757:            WRITE(MYUNIT,'(A,I6,A,E20.10)') &
758:                 'finalio> Min ',J1,' has total mass ', EDUMMY758:                 'finalio> Min ',J1,' has total mass ', EDUMMY
759:         ENDIF759:         ENDIF
760:         !760:         !
761:         IF(STRESST) CALL CALC_STRESS(QMINP(J1,1:3*NATOMS),EDUMMY)761:         IF(STRESST) CALL CALC_STRESS(QMINP(J1,1:3*NATOMS),EDUMMY)
762:         !762:         !
763:         DO J2=1,NSPECIES(0) !ds656> Should not exceed 10763:         DO J2=1,NSPECIES(0) !ds656> Should not exceed 10
764:            IF(SPECLABELST) THEN764:            IF(SPECLABELST) THEN
842:   !842:   !
843:   ! End of loop over dump to file lowest or equivalent. 843:   ! End of loop over dump to file lowest or equivalent. 
844:   !844:   !
845:   CLOSE(MYUNIT2)845:   CLOSE(MYUNIT2)
846:   IF (CSMT.AND.(.NOT.SYMMETRIZECSM)) CLOSE(MYUNIT3)846:   IF (CSMT.AND.(.NOT.SYMMETRIZECSM)) CLOSE(MYUNIT3)
847:   !847:   !
848:   !     csw34> New loop for dumping interaction energy files if A9INTE is specified848:   !     csw34> New loop for dumping interaction energy files if A9INTE is specified
849:   !     Added the missing IF block to test for A9INTE 9/12/09 DJW849:   !     Added the missing IF block to test for A9INTE 9/12/09 DJW
850:   !850:   !
851:   IF (A9INTET) THEN851:   IF (A9INTET) THEN
852:      IF (MPIT.OR.GBHT) THEN852:      IF (MPIT) THEN
853:         WRITE (ISTR, '(I10)') MYUNIT-22980+1853:         WRITE (ISTR, '(I10)') MYUNIT-22980+1
854:         MYUNIT2=(MYUNIT-22980+1)+100854:         MYUNIT2=(MYUNIT-22980+1)+100
855:         MYFILENAME2="intelowest."//trim(adjustl(istr))855:         MYFILENAME2="intelowest."//trim(adjustl(istr))
856:         OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")856:         OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")
857:      ELSE857:      ELSE
858:         MYUNIT2=25858:         MYUNIT2=25
859:         OPEN(MYUNIT2,FILE='intelowest',STATUS='UNKNOWN')859:         OPEN(MYUNIT2,FILE='intelowest',STATUS='UNKNOWN')
860:      ENDIF860:      ENDIF
861:   ENDIF861:   ENDIF
862:   !862:   !
954:                 SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))*COS(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+2)), &954:                 SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))*COS(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+2)), &
955:                 SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))*SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+2)), &955:                 SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))*SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+2)), &
956:                 COS(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))956:                 COS(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))
957:         ENDDO957:         ENDDO
958:      ENDDO958:      ENDDO
959:      CLOSE(LUNIT)959:      CLOSE(LUNIT)
960:      960:      
961:   ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT.OR.PYT) THEN961:   ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT.OR.PYT) THEN
962:      DO J1=1,NSAVE962:      DO J1=1,NSAVE
963:         WRITE(J1CHAR2,'(I3)') J1963:         WRITE(J1CHAR2,'(I3)') J1
964:         IF (MPIT.OR.GBHT) THEN964:         IF (MPIT) THEN
965:            WRITE (ISTR, '(I10)') MYNODE+1965:            WRITE (ISTR, '(I10)') MYNODE+1
966:            MYUNIT2=GETUNIT()966:            MYUNIT2=GETUNIT()
967:            MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))//'.'//trim(adjustl(istr))967:            MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))//'.'//trim(adjustl(istr))
968:            OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")968:            OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")
969:         ELSE969:         ELSE
970:            MYUNIT2=GETUNIT()970:            MYUNIT2=GETUNIT()
971:            MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))971:            MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))
972:            OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")972:            OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")
973:         END IF973:         END IF
974:         974:         
977:         ENDDO977:         ENDDO
978:         978:         
979:         CLOSE(MYUNIT2)979:         CLOSE(MYUNIT2)
980:         980:         
981:      ENDDO981:      ENDDO
982:      !982:      !
983:      !  Write out lowest NSAVE structures to xmakemol ellipsoid format for983:      !  Write out lowest NSAVE structures to xmakemol ellipsoid format for
984:      !  clusters of ellipsoids of revolution984:      !  clusters of ellipsoids of revolution
985:      !985:      !
986:      986:      
987:      IF (MPIT.OR.GBHT) THEN987:      IF (MPIT) THEN
988:         WRITE (ISTR, '(I10)') MYNODE+1988:         WRITE (ISTR, '(I10)') MYNODE+1
989:         MYUNIT2=GETUNIT()989:         MYUNIT2=GETUNIT()
990:         MYFILENAME2="ellipsoid."//trim(adjustl(istr))//".xyz"990:         MYFILENAME2="ellipsoid."//trim(adjustl(istr))//".xyz"
991:         OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")991:         OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")
992:      ELSE992:      ELSE
993:         MYUNIT2=GETUNIT()993:         MYUNIT2=GETUNIT()
994:         OPEN(MYUNIT2,FILE='ellipsoid.xyz',STATUS='UNKNOWN')994:         OPEN(MYUNIT2,FILE='ellipsoid.xyz',STATUS='UNKNOWN')
995:      ENDIF995:      ENDIF
996:      996:      
997:      IF(PYT) THEN997:      IF(PYT) THEN


r31263/initialization.f90 2016-10-12 11:30:25.624806187 +0100 r31262/initialization.f90 2016-10-12 11:30:31.524885391 +0100
 22: !   that at this point the keywords have been read and the coords have been 22: !   that at this point the keywords have been read and the coords have been
 23: !   loaded. 23: !   loaded.
 24:  24: 
 25: SUBROUTINE INITIALIZATIONS() 25: SUBROUTINE INITIALIZATIONS()
 26: USE COMMONS 26: USE COMMONS
 27: USE QMODULE 27: USE QMODULE
 28: USE CLASS_OVERLAP 28: USE CLASS_OVERLAP
 29: USE BGUPMOD ! In BGupta.f90 29: USE BGUPMOD ! In BGupta.f90
 30: USE POT_PARAMS  30: USE POT_PARAMS 
 31: USE POLIRMOD, ONLY: POLIRINIT 31: USE POLIRMOD, ONLY: POLIRINIT
 32: !USE HOMOREFMOD ! in homoref_addons.f90 32: USE HOMOREFMOD ! in homoref_addons.f90
 33: IMPLICIT NONE 33: IMPLICIT NONE
 34: INTEGER LUNIT, GETUNIT 34: INTEGER LUNIT, GETUNIT
 35: integer j1, j2, j6 35: integer j1, j2, j6
 36: DOUBLE PRECISION TEMPCOORDS(3*NATOMS) 36: DOUBLE PRECISION TEMPCOORDS(3*NATOMS)
 37: DOUBLE PRECISION GRAD(NATOMS*3), ENERGY, DUMMY 37: DOUBLE PRECISION GRAD(NATOMS*3), ENERGY, DUMMY
 38: LOGICAL GTEST 38: LOGICAL GTEST
 39:  39: 
 40: IF (ONEDAPBCT.OR.ONEDPBCT) THEN 40: IF (ONEDAPBCT.OR.ONEDPBCT) THEN
 41:    ALLOCATE(XYPHI(3*NATOMS)) 41:    ALLOCATE(XYPHI(3*NATOMS))
 42:    LUNIT=GETUNIT() 42:    LUNIT=GETUNIT()
294:          !294:          !
295:          DUMMY = DUMMY/CUTOFF295:          DUMMY = DUMMY/CUTOFF
296:          CUTC_ATT(J1,J2) = -DBLE(MSC_M(J1,J2)*(MSC_M(J1,J2)+1))*DUMMY296:          CUTC_ATT(J1,J2) = -DBLE(MSC_M(J1,J2)*(MSC_M(J1,J2)+1))*DUMMY
297:          CUTC_ATT(J2,J1) = CUTC_ATT(J1,J2)297:          CUTC_ATT(J2,J1) = CUTC_ATT(J1,J2)
298:          !298:          !
299:       ENDDO299:       ENDDO
300:    ENDDO300:    ENDDO
301:    !301:    !
302: ENDIF302: ENDIF
303: !303: !
 304: IF(HOMOREFT) THEN
 305:    ALLOCATE(NNBRS(NATOMS,NSPECIES(0),0:NNMAX))
 306:    ALLOCATE(ANNHIST(NSPECIES(0),NSPECIES(0),-1:NNMAX))
 307:    ALLOCATE(NNHIST(NSPECIES(0),NSPECIES(0),-1:NNMAX))
 308:    ALLOCATE(ANNHIST_MEAN(NSPECIES(0),NSPECIES(0)))
 309:    ALLOCATE(ANNHIST_VAR(NSPECIES(0),NSPECIES(0)))
 310:    ALLOCATE(NNBOND_WEIGHTS(NSPECIES(0),NSPECIES(0)))
 311:    ALLOCATE(IFLIPE(1:NATOMS))
 312: ENDIF
 313: !
304: ! ds656 > Allocate book-keeping arrays for listing314: ! ds656 > Allocate book-keeping arrays for listing
305: ! atoms and their neighbours.  315: ! atoms and their neighbours.  
306: !316: !
307: ALLOCATE(ATOMLISTS(NSPECIES(0),3,0:NATOMS))317: ALLOCATE(ATOMLISTS(NSPECIES(0),3,0:NATOMS))
308: ALLOCATE(INVATOMLISTS(NATOMS,3))318: ALLOCATE(INVATOMLISTS(NATOMS,3))
309: !319: !
310: IF(QALCS_SURFT .OR. QALCS_SYMT) THEN320: IF(QALCS_SURFT .OR. QALCS_SYMT) THEN
311:    NNLIST_MAX = 12321:    NNLIST_MAX = 12
312:    ALLOCATE(NNLISTS(1:NATOMS,0:NNLIST_MAX))322:    ALLOCATE(NNLISTS(1:NATOMS,0:NNLIST_MAX))
313:    IF(QALCS_SURFT) THEN323:    IF(QALCS_SURFT) THEN


r31263/io1.f 2016-10-12 11:30:25.880809635 +0100 r31262/io1.f 2016-10-12 11:30:31.784888879 +0100
392:             ENDIF392:             ENDIF
393:             DO J1=NTYPEA+1,NATOMS393:             DO J1=NTYPEA+1,NATOMS
394:                IF (ABS(ATMASS(J1)-ATMASS(NTYPEA+1)).GT.1.0D-10) THEN394:                IF (ABS(ATMASS(J1)-ATMASS(NTYPEA+1)).GT.1.0D-10) THEN
395:                   WRITE(MYUNIT,'(A,I6,A,2G20.10)') 'io1> WARNING *** atom masses for NTYPEA+1 and ',J1,' are ',395:                   WRITE(MYUNIT,'(A,I6,A,2G20.10)') 'io1> WARNING *** atom masses for NTYPEA+1 and ',J1,' are ',
396:      &                 ATMASS(NTYPEA+1),ATMASS(J1)396:      &                 ATMASS(NTYPEA+1),ATMASS(J1)
397:                ENDIF397:                ENDIF
398:             ENDDO398:             ENDDO
399:          ELSE399:          ELSE
400:             WRITE(MYUNIT,'(A)') 'Atom masses will be set to species indices.'400:             WRITE(MYUNIT,'(A)') 'Atom masses will be set to species indices.'
401:          ENDIF !...<ds656401:          ENDIF !...<ds656
 402:       ELSE IF (BASWAP) THEN
 403:          WRITE(MYUNIT,'(A,I6)') 'Atom identity swaps will be performed. NWAIT= ', BASWAP_NWAIT
 404:          WRITE(MYUNIT,'(A,F10.6)') 'Maximum fraction of swap moves between A and
 405:      1     B type particles:  ', BASWAP_FRAC
402:       ELSE IF (DZTEST) THEN406:       ELSE IF (DZTEST) THEN
403:          WRITE(MYUNIT,'(I4,A,F10.4)') NATOMS,' Dzugutov atoms'407:          WRITE(MYUNIT,'(I4,A,F10.4)') NATOMS,' Dzugutov atoms'
404:       ELSE IF (ZETT1) THEN408:       ELSE IF (ZETT1) THEN
405:          WRITE(MYUNIT,'(I4,A,F10.4)') NATOMS,' Zetterling type 1 atoms'409:          WRITE(MYUNIT,'(I4,A,F10.4)') NATOMS,' Zetterling type 1 atoms'
406:       ELSE IF (ZETT2) THEN410:       ELSE IF (ZETT2) THEN
407:          WRITE(MYUNIT,'(I4,A,F10.4)') NATOMS,' Zetterling type 2 atoms'411:          WRITE(MYUNIT,'(I4,A,F10.4)') NATOMS,' Zetterling type 2 atoms'
408:       ELSE IF (PACHECO) THEN412:       ELSE IF (PACHECO) THEN
409:          WRITE(MYUNIT,'(I4,A,F10.4)') NATOMS,' Pacheco-Ramalho C60 molecules'413:          WRITE(MYUNIT,'(I4,A,F10.4)') NATOMS,' Pacheco-Ramalho C60 molecules'
410:       ELSE IF (MODEL1T) THEN414:       ELSE IF (MODEL1T) THEN
411:          WRITE(MYUNIT,'(A)') ' 1D landscape model 1'415:          WRITE(MYUNIT,'(A)') ' 1D landscape model 1'


r31263/keywords.f 2016-10-12 11:30:26.140813124 +0100 r31262/keywords.f 2016-10-12 11:30:32.064892637 +0100
488:       LSWAPS_TEMP = 0.1488:       LSWAPS_TEMP = 0.1
489:       LSWAPS_TFACTOR = 1.0489:       LSWAPS_TFACTOR = 1.0
490:       LSWAPS_NUP = 1490:       LSWAPS_NUP = 1
491:       LFLIPST=.FALSE.491:       LFLIPST=.FALSE.
492:       LFLIPS_RESET=.FALSE.492:       LFLIPS_RESET=.FALSE.
493:       LFLIPS_N1 = 1493:       LFLIPS_N1 = 1
494:       LFLIPS_N2 = 1494:       LFLIPS_N2 = 1
495:       LFLIPS_TEMP = 0.1495:       LFLIPS_TEMP = 0.1
496:       LFLIPS_TFACTOR = 1.0496:       LFLIPS_TFACTOR = 1.0
497:       LFLIPS_NUP = 1497:       LFLIPS_NUP = 1
 498:       BASWAP=.FALSE.
 499:       BASWAP_FRAC=0.0D0
 500:       !BASWAP_TEMP=TEMP(1) ! no longer used
 501:       BASWAP_NWAIT=0
 502:       BASWAP_NMAX=0
 503:       BASWAPTEST=.FALSE.
498:       NTYPEA = NATOMS ! ds656> before this was initialised to 0504:       NTYPEA = NATOMS ! ds656> before this was initialised to 0
499:       BGUPTANAME1 = 'Au'505:       BGUPTANAME1 = 'Au'
500:       BGUPTANAME2 = 'Ag'506:       BGUPTANAME2 = 'Ag'
 507:       HOMOREFT = .FALSE.
 508:       HOMOREFTEST = .FALSE.
 509:       HOMOREFCHECK = .FALSE.
 510:       HOMOREF_LSMODE = 0
 511:       HOMOREF_FGMODE = 0
 512:       HOMOREF_NCYCLES = 0
 513:       HOMOREF_NFMAX = 3
 514:       HOMOREF_NSMAX = 1
 515:       HOMOREF_AUXT = .FALSE.
 516:       HOMOREF_AUX_NSWAPS = 0
 517:       HOMOREF_AUX_TEMP = 1.0D0
 518:       HOMOREF_AUX_FACTOR = 1.0D0
 519:       HOMOREF_AUX_NNCUT = 1.5D0
 520:       HOMOREF_BHT = .FALSE.
 521:       HOMOREF_BH_NSWAPMAX = 0
 522:       HOMOREF_BH_NDUDMAX = 0
 523:       HOMOREF_BH_TEMP = 1.0D0
 524:       HOMOREF_BH_FACTOR = 1.0D0
501:       QALCST = .FALSE.525:       QALCST = .FALSE.
502:       QALCSV = .FALSE.526:       QALCSV = .FALSE.
503:       QALCSMODE = 0527:       QALCSMODE = 0
504:       QALCS_PARAM = 0528:       QALCS_PARAM = 0
505:       QALCS_SURFT = .FALSE.529:       QALCS_SURFT = .FALSE.
506:       QALCS_SYMT = .FALSE.530:       QALCS_SYMT = .FALSE.
507:       QALCS_SYM_MINCORESIZE = 5531:       QALCS_SYM_MINCORESIZE = 5
508:       SPECLABELST = .FALSE.532:       SPECLABELST = .FALSE.
509:       RANDMULTIPERMT = .FALSE.533:       RANDMULTIPERMT = .FALSE.
510:       RANDMULTIPERM_STEP=1534:       RANDMULTIPERM_STEP=1
1081: ! jdf43> MBPOL1105: ! jdf43> MBPOL
1082:       MBPOLT=.FALSE.1106:       MBPOLT=.FALSE.
1083:       MOLECULART=.FALSE.1107:       MOLECULART=.FALSE.
1084: ! jdf43>1108: ! jdf43>
1085:       REPMATCHT=.FALSE.1109:       REPMATCHT=.FALSE.
1086:       1110:       
1087:       UNIFORMMOVE=.FALSE.1111:       UNIFORMMOVE=.FALSE.
1088:       ORBITTOL=1.0D-31112:       ORBITTOL=1.0D-3
1089:       NOINVERSION=.FALSE.1113:       NOINVERSION=.FALSE.
1090: !1114: !
1091: ! ds656> Parallelised generalised basin-hopping  
1092:       GBHT=.FALSE. 
1093:  
1094: ! General mixed LJ systems1115: ! General mixed LJ systems
1095:       GLJT=.FALSE.1116:       GLJT=.FALSE.
1096: ! ds656> Multicomponent LJ system (different implementation to GLJ!)1117: ! ds656> Multicomponent LJ system (different implementation to GLJ!)
1097:       MLJT=.FALSE.1118:       MLJT=.FALSE.
1098: 1119: 
1099: ! khs26> Free energy basin-hopping stuff1120: ! khs26> Free energy basin-hopping stuff
1100:       FEBHT = .FALSE.1121:       FEBHT = .FALSE.
1101:       FETEMP = 0.0D01122:       FETEMP = 0.0D0
1102: ! khs26> This requires a minimum separation between the zero and non-zero1123: ! khs26> This requires a minimum separation between the zero and non-zero
1103: ! eigenvalues for runs using free energy basin-hopping. This is based on 1124: ! eigenvalues for runs using free energy basin-hopping. This is based on 
1441:          GLJEPS(1,2) = EPSAB1462:          GLJEPS(1,2) = EPSAB
1442:          GLJEPS(2,1) = EPSAB ! impose symmetry1463:          GLJEPS(2,1) = EPSAB ! impose symmetry
1443:          GLJSIG(1,1) = 1.0D01464:          GLJSIG(1,1) = 1.0D0
1444:          GLJSIG(2,2) = SIGBB1465:          GLJSIG(2,2) = SIGBB
1445:          GLJSIG(1,2) = SIGAB1466:          GLJSIG(1,2) = SIGAB
1446:          GLJSIG(2,1) = SIGAB ! impose symmetry1467:          GLJSIG(2,1) = SIGAB ! impose symmetry
1447: 1468: 
1448:       ELSE IF (WORD .EQ. 'PHI4MODEL') THEN1469:       ELSE IF (WORD .EQ. 'PHI4MODEL') THEN
1449:          PHI4MODELT = .TRUE.1470:          PHI4MODELT = .TRUE.
1450: 1471: 
1451:       ELSE IF (WORD .EQ. 'GBH') THEN1472: ! ds656> Homotop refinement for binary systems
1452:          GBHT=.TRUE.1473:       ELSE IF (WORD .EQ. 'HOMOREF') THEN
1453: 1474:          HOMOREFT = .TRUE.
 1475:          IF(NITEMS >= 5) THEN
 1476:             CALL READI(HOMOREF_LSMODE)
 1477:             CALL READI(HOMOREF_FGMODE)
 1478:             CALL READI(HOMOREF_NCYCLES)
 1479:             CALL READI(HOMOREF_NFMAX)
 1480:             IF(NITEMS > 5) CALL READI(HOMOREF_NSMAX)
 1481:             WRITE(MYUNIT,'(A,I2,A,I2,A,I5,A,I2,A,I2)') 
 1482:      2           'keywords> HOMOREF with LSMODE=', HOMOREF_LSMODE, 
 1483:      3           ', FGMODE=', HOMOREF_FGMODE, 
 1484:      4           ', NCYCLES=', HOMOREF_NCYCLES,
 1485:      5           ', NFMAX=', HOMOREF_NFMAX,
 1486:      6           ' and NSMAX=', HOMOREF_NSMAX
 1487:          ELSE
 1488:             WRITE(MYUNIT,'(A)') 
 1489:      2           "keywords> HOMOREF requires at least 5 arguments!"
 1490:             STOP
 1491:          ENDIF
 1492:       ELSE IF (WORD .EQ. 'HOMOREF_AUX') THEN
 1493:          HOMOREF_AUXT = .TRUE.
 1494:          IF(NITEMS .EQ. 5) THEN
 1495:             CALL READI(HOMOREF_AUX_NSWAPS)
 1496:             CALL READF(HOMOREF_AUX_TEMP)
 1497:             CALL READF(HOMOREF_AUX_FACTOR)
 1498:             CALL READF(HOMOREF_AUX_NNCUT)
 1499:             WRITE(MYUNIT,'(A,I4,A,F6.3,A,F6.3,A,F8.4)') 
 1500:      2           'keywords> HOMOREF_AUX with NSWAPS=',
 1501:      3           HOMOREF_AUX_NSWAPS,', TEMP=', 
 1502:      4           HOMOREF_AUX_TEMP,', FACTOR=',
 1503:      5           HOMOREF_AUX_FACTOR,' and NNCUT=',
 1504:      6           HOMOREF_AUX_NNCUT
 1505:          ELSE
 1506:             WRITE(MYUNIT,'(A)') "keywords> HOMOREF_AUX takes 4 arguments!"
 1507:             STOP            
 1508:          ENDIF
 1509:       ELSE IF (WORD .EQ. 'HOMOREF_BH') THEN
 1510:          HOMOREF_BHT = .TRUE.
 1511:          IF(NITEMS .EQ. 5) THEN
 1512:             CALL READI(HOMOREF_BH_NSWAPMAX)
 1513:             CALL READI(HOMOREF_BH_NDUDMAX)
 1514:             CALL READF(HOMOREF_BH_TEMP)
 1515:             CALL READF(HOMOREF_BH_FACTOR)
 1516:          ELSE
 1517:             WRITE(MYUNIT,'(A)') "keywords> HOMOREF_BH takes 4 arguments!"
 1518:             STOP            
 1519:          ENDIF
 1520:       ELSE IF (WORD .EQ. 'HOMOREFTEST') THEN
 1521:          HOMOREFTEST = .TRUE.
 1522:       ELSE IF (WORD .EQ. 'HOMOREFCHECK') THEN
 1523:          HOMOREFCHECK = .TRUE.
1454: ! ds656> Quanech-Assisted Local Combinatorial Search1524: ! ds656> Quanech-Assisted Local Combinatorial Search
1455:       ELSE IF (WORD .EQ. 'QALCS') THEN1525:       ELSE IF (WORD .EQ. 'QALCS') THEN
1456:          QALCST = .TRUE.1526:          QALCST = .TRUE.
1457:          IF(NITEMS > 1) THEN1527:          IF(NITEMS > 1) THEN
1458:             CALL READI(QALCSMODE)1528:             CALL READI(QALCSMODE)
1459:             IF(NITEMS > 2) THEN1529:             IF(NITEMS > 2) THEN
1460:                CALL READI(QALCS_PARAM)1530:                CALL READI(QALCS_PARAM)
1461:             ENDIF1531:             ENDIF
1462:          ELSE1532:          ELSE
1463:             WRITE(MYUNIT,'(A)') 'keywords> QALCS mode unspecified!'1533:             WRITE(MYUNIT,'(A)') 'keywords> QALCS mode unspecified!'
1524:          CALL READF(PAB)1594:          CALL READF(PAB)
1525:          CALL READF(QAB)1595:          CALL READF(QAB)
1526:          CALL READF(ZAB)1596:          CALL READF(ZAB)
1527:          CALL READF(R0AB)1597:          CALL READF(R0AB)
1528:        ELSE IF (WORD.EQ.'BGUPTATBB') THEN1598:        ELSE IF (WORD.EQ.'BGUPTATBB') THEN
1529:          CALL READF(ABB)1599:          CALL READF(ABB)
1530:          CALL READF(PBB)1600:          CALL READF(PBB)
1531:          CALL READF(QBB)1601:          CALL READF(QBB)
1532:          CALL READF(ZBB)1602:          CALL READF(ZBB)
1533:          CALL READF(R0BB)1603:          CALL READF(R0BB)
 1604:        ELSE IF (WORD.EQ.'BASWAPTEST') THEN
 1605:           BASWAPTEST=.TRUE.
 1606:        ELSE IF (WORD.EQ.'FIXLABELS') THEN
 1607:           FIXLABELS=.TRUE.
1534:        ELSE IF (WORD.EQ.'LSWAPS') THEN1608:        ELSE IF (WORD.EQ.'LSWAPS') THEN
1535:          LSWAPST=.TRUE.1609:          LSWAPST=.TRUE.
1536:          CALL READI(LSWAPS_N1)1610:          CALL READI(LSWAPS_N1)
1537:          CALL READI(LSWAPS_N2)1611:          CALL READI(LSWAPS_N2)
1538:          CALL READF(LSWAPS_TEMP)1612:          CALL READF(LSWAPS_TEMP)
1539:          IF(NITEMS > 4) CALL READF(LSWAPS_TFACTOR)1613:          IF(NITEMS > 4) CALL READF(LSWAPS_TFACTOR)
1540:          IF(NITEMS > 5) CALL READI(LSWAPS_NUP)1614:          IF(NITEMS > 5) CALL READI(LSWAPS_NUP)
1541:       ELSE IF (WORD.EQ.'LFLIPS') THEN1615:       ELSE IF (WORD.EQ.'LFLIPS') THEN
1542:          LFLIPST=.TRUE.1616:          LFLIPST=.TRUE.
1543:          CALL READI(LFLIPS_N1)1617:          CALL READI(LFLIPS_N1)
1544:          CALL READI(LFLIPS_N2)1618:          CALL READI(LFLIPS_N2)
1545:          CALL READF(LFLIPS_TEMP)1619:          CALL READF(LFLIPS_TEMP)
1546:          IF(NITEMS > 4) CALL READF(LFLIPS_TFACTOR)1620:          IF(NITEMS > 4) CALL READF(LFLIPS_TFACTOR)
1547:          IF(NITEMS > 5) CALL READI(LFLIPS_NUP)1621:          IF(NITEMS > 5) CALL READI(LFLIPS_NUP)
1548:       ELSE IF (WORD.EQ.'LFLIPS_RESET') THEN1622:       ELSE IF (WORD.EQ.'LFLIPS_RESET') THEN
1549:          LFLIPS_RESET=.TRUE.1623:          LFLIPS_RESET=.TRUE.
 1624:       ELSE IF (WORD.EQ.'BASWAP') THEN
 1625:          BASWAP=.TRUE.
 1626:          CALL READI(BASWAP_NWAIT)
 1627:          !CALL READF(BASWAP_TEMP) ! no longer used!
 1628:          CALL READF(BASWAP_FRAC)
 1629:          CALL READI(BASWAP_NMAX)
1550:       ELSE IF (WORD.EQ.'RANDPERM') THEN1630:       ELSE IF (WORD.EQ.'RANDPERM') THEN
1551:          RANDPERMT=.TRUE.1631:          RANDPERMT=.TRUE.
1552:          IF(NITEMS > 1) THEN1632:          IF(NITEMS > 1) THEN
1553:             CALL READI(RANDPERM_STEP)1633:             CALL READI(RANDPERM_STEP)
1554:          ENDIF1634:          ENDIF
1555:       ELSE IF (WORD.EQ.'RANDMULTIPERM') THEN1635:       ELSE IF (WORD.EQ.'RANDMULTIPERM') THEN
1556:          RANDMULTIPERMT=.TRUE.1636:          RANDMULTIPERMT=.TRUE.
1557:          IF(NITEMS > 1) THEN1637:          IF(NITEMS > 1) THEN
1558:             CALL READI(RANDMULTIPERM_STEP)1638:             CALL READI(RANDMULTIPERM_STEP)
1559:          ENDIF1639:          ENDIF


r31263/main.F 2016-10-12 11:30:26.400816606 +0100 r31262/main.F 2016-10-12 11:30:32.324896128 +0100
 23:       USE COMMONS 23:       USE COMMONS
 24:       USE QMODULE 24:       USE QMODULE
 25:       USE PERMU 25:       USE PERMU
 26:       USE F1COM 26:       USE F1COM
 27:       USE MODAMBER 27:       USE MODAMBER
 28:       USE MODAMBER9, only : AMBFINALIO_NODE,MDCRD_UNIT,MDINFO_UNIT,AMBPDB_UNIT, ATMASS1 28:       USE MODAMBER9, only : AMBFINALIO_NODE,MDCRD_UNIT,MDINFO_UNIT,AMBPDB_UNIT, ATMASS1
 29:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_MASSES 29:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_MASSES
 30:       USE MODCHARMM 30:       USE MODCHARMM
 31:       USE PORFUNCS 31:       USE PORFUNCS
 32:       USE TWIST_MOD 32:       USE TWIST_MOD
 33:       !USE HOMOREFMOD 33:       USE HOMOREFMOD
 34:       USE GENRIGID, only : RIGIDINIT, GENRIGID_READ_FROM_FILE, DEGFREEDOMS 34:       USE GENRIGID, only : RIGIDINIT, GENRIGID_READ_FROM_FILE, DEGFREEDOMS
 35:       USE MULTIPOT, only: MULTIPOT_INITIALISE 35:       USE MULTIPOT, only: MULTIPOT_INITIALISE
 36:       USE GAUSS_MOD, ONLY: KEGEN 36:       USE GAUSS_MOD, ONLY: KEGEN
 37:       USE MOLECULAR_DYNAMICS, ONLY : MDT, MD_NSTEPS, MD_NWAIT, MDRUN 37:       USE MOLECULAR_DYNAMICS, ONLY : MDT, MD_NSTEPS, MD_NWAIT, MDRUN
 38:  38: 
 39:       IMPLICIT NONE 39:       IMPLICIT NONE
 40:       !EXTERNAL READ_CMD_ARGS 40:       !EXTERNAL READ_CMD_ARGS
 41: #ifdef MPI 41: #ifdef MPI
 42:       INCLUDE 'mpif.h' 42:       INCLUDE 'mpif.h'
 43: #endif 43: #endif
 75:       IF (TRIM(ADJUSTL(INFILE)).EQ.'') INFILE='output' 75:       IF (TRIM(ADJUSTL(INFILE)).EQ.'') INFILE='output'
 76:       IF (NPAR.GT.1) THEN  76:       IF (NPAR.GT.1) THEN 
 77:          MYFILENAME=TRIM(ADJUSTL(INFILE))//"."//TRIM(ADJUSTL(ISTR)) 77:          MYFILENAME=TRIM(ADJUSTL(INFILE))//"."//TRIM(ADJUSTL(ISTR))
 78:       ELSE 78:       ELSE
 79:          MYFILENAME=TRIM(ADJUSTL(INFILE)) 79:          MYFILENAME=TRIM(ADJUSTL(INFILE))
 80:       ENDIF 80:       ENDIF
 81:       OPEN(MYUNIT,FILE=MYFILENAME, STATUS="unknown", form="formatted") 81:       OPEN(MYUNIT,FILE=MYFILENAME, STATUS="unknown", form="formatted")
 82:       WRITE(MYUNIT, '(A,I10,A,I10)') "Starting parallel execution: Processor", mynode+1, " of ",NPAR 82:       WRITE(MYUNIT, '(A,I10,A,I10)') "Starting parallel execution: Processor", mynode+1, " of ",NPAR
 83: #else 83: #else
 84:       NPAR=1 84:       NPAR=1
 85:       NPAR_GBH=1 
 86:       MYNODE=0 85:       MYNODE=0
 87:       MYUNIT=22979+1 86:       MYUNIT=22979+1
 88:       MDCRD_UNIT=20000 87:       MDCRD_UNIT=20000
 89:       MDINFO_UNIT=21000 88:       MDINFO_UNIT=21000
 90:       AMBPDB_UNIT=18000 89:       AMBPDB_UNIT=18000
 91:       AMBFINALIO_NODE=1 90:       AMBFINALIO_NODE=1
 92:       IF (TRIM(ADJUSTL(INFILE)).EQ.'') THEN 91:       IF (TRIM(ADJUSTL(INFILE)).EQ.'') THEN
 93:          MYFILENAME="output" 92:          MYFILENAME="output"
 94:       ELSE 93:       ELSE
 95:          MYFILENAME=TRIM(ADJUSTL(INFILE)) 94:          MYFILENAME=TRIM(ADJUSTL(INFILE))
110: !         WRITE(MYUNIT, '(A)') "NOTE: Binary only licensed for use at AlgoSB2015"109: !         WRITE(MYUNIT, '(A)') "NOTE: Binary only licensed for use at AlgoSB2015"
111: !      ELSE110: !      ELSE
112: ! if after the workshop has ended, stop execution with error111: ! if after the workshop has ended, stop execution with error
113: !         WRITE(MYUNIT, '(A)') "ERROR: license missing, contact dw34@cam.ac.uk. Stopping."112: !         WRITE(MYUNIT, '(A)') "ERROR: license missing, contact dw34@cam.ac.uk. Stopping."
114: !         STOP   113: !         STOP   
115: !      END IF114: !      END IF
116: 115: 
117: ! Add the GMIN version to the output - helps bug hunting :)116: ! Add the GMIN version to the output - helps bug hunting :)
118:       !WRITE(MYUNIT, '(A,I5)') 'GMIN version r',VERSIONTEMP117:       !WRITE(MYUNIT, '(A,I5)') 'GMIN version r',VERSIONTEMP
119: !     CALL DISPLAY_VERSION(MYUNIT)118: !     CALL DISPLAY_VERSION(MYUNIT)
120:       CALL COUNTATOMS(MYUNIT, NPAR, GCBHT, GCMU, GCNATOMS, GBHT)119:       CALL COUNTATOMS(MYUNIT, NPAR, GCBHT, GCMU, GCNATOMS)
121:       IF(GBHT) THEN 
122:          NPAR_GBH=NPAR 
123:          NPAR=1 ! ds656> To prevent the usual memory allocation. 
124:       ENDIF 
125:       IF (CHRMMT) WRITE(MYUNIT,'(A,I8)') 'main> MAXAIM parameter for CHARMM MXATMS=',MXATMS120:       IF (CHRMMT) WRITE(MYUNIT,'(A,I8)') 'main> MAXAIM parameter for CHARMM MXATMS=',MXATMS
126: !121: !
127: ! NATOMS is set to NUMBER_OF_ATOMS in modcommonunit122: ! NATOMS is set to NUMBER_OF_ATOMS in modcommonunit
128: !123: !
129:       IF (MLPB3T) THEN124:       IF (MLPB3T) THEN
130:          IF (NUMBER_OF_ATOMS.NE.MLPHIDDEN*(MLPIN+MLPOUT)+1) THEN125:          IF (NUMBER_OF_ATOMS.NE.MLPHIDDEN*(MLPIN+MLPOUT)+1) THEN
131:             WRITE(MYUNIT,'(A,I8,A,I8)') 'main> ERROR *** number of coordinates in coords is ',NUMBER_OF_ATOMS,126:             WRITE(MYUNIT,'(A,I8,A,I8)') 'main> ERROR *** number of coordinates in coords is ',NUMBER_OF_ATOMS,
132:      &                         ' should be ',MLPHIDDEN*(MLPIN+MLPOUT) +1127:      &                         ' should be ',MLPHIDDEN*(MLPIN+MLPOUT) +1
133:             STOP128:             STOP
134:          ENDIF129:          ENDIF
684:       IF (ALLOCATED(ATOMLISTS)) DEALLOCATE(ATOMLISTS)679:       IF (ALLOCATED(ATOMLISTS)) DEALLOCATE(ATOMLISTS)
685:       IF (ALLOCATED(INVATOMLISTS)) DEALLOCATE(INVATOMLISTS)680:       IF (ALLOCATED(INVATOMLISTS)) DEALLOCATE(INVATOMLISTS)
686:       IF (ALLOCATED(NBRLISTS)) DEALLOCATE(NBRLISTS)681:       IF (ALLOCATED(NBRLISTS)) DEALLOCATE(NBRLISTS)
687:       IF (ALLOCATED(NNLISTS)) DEALLOCATE(NNLISTS)682:       IF (ALLOCATED(NNLISTS)) DEALLOCATE(NNLISTS)
688:       IF (ALLOCATED(VSITES)) DEALLOCATE(VSITES)683:       IF (ALLOCATED(VSITES)) DEALLOCATE(VSITES)
689:       IF (ALLOCATED(MIEF_SIG)) DEALLOCATE(MIEF_SIG)684:       IF (ALLOCATED(MIEF_SIG)) DEALLOCATE(MIEF_SIG)
690:       IF (ALLOCATED(MIEF_EPS)) DEALLOCATE(MIEF_EPS)685:       IF (ALLOCATED(MIEF_EPS)) DEALLOCATE(MIEF_EPS)
691:       IF (ALLOCATED(MIEF_SITES)) DEALLOCATE(MIEF_SITES)686:       IF (ALLOCATED(MIEF_SITES)) DEALLOCATE(MIEF_SITES)
692:       IF (ALLOCATED(MIEF_U_RCUT)) DEALLOCATE(MIEF_U_RCUT)687:       IF (ALLOCATED(MIEF_U_RCUT)) DEALLOCATE(MIEF_U_RCUT)
693:       IF (ALLOCATED(MIEF_DUDR_RCUT)) DEALLOCATE(MIEF_DUDR_RCUT)688:       IF (ALLOCATED(MIEF_DUDR_RCUT)) DEALLOCATE(MIEF_DUDR_RCUT)
 689:       IF (ALLOCATED(NNBRS)) DEALLOCATE(NNBRS)
 690:       IF (ALLOCATED(NNHIST)) DEALLOCATE(NNHIST)
 691:       IF (ALLOCATED(ANNHIST)) DEALLOCATE(ANNHIST)
 692:       IF (ALLOCATED(ANNHIST_MEAN)) DEALLOCATE(ANNHIST_MEAN)
 693:       IF (ALLOCATED(ANNHIST_VAR)) DEALLOCATE(ANNHIST_VAR)
 694:       IF (ALLOCATED(NNBOND_WEIGHTS)) DEALLOCATE(NNBOND_WEIGHTS)
 695:       IF (ALLOCATED(IFLIPE)) DEALLOCATE(IFLIPE)
694:       IF (ALLOCATED(OBJ)) DEALLOCATE(OBJ)696:       IF (ALLOCATED(OBJ)) DEALLOCATE(OBJ)
695:       IF (ALLOCATED(TWIST_K)) DEALLOCATE(TWIST_K)697:       IF (ALLOCATED(TWIST_K)) DEALLOCATE(TWIST_K)
696:       IF (ALLOCATED(TWIST_THETA0)) DEALLOCATE(TWIST_THETA0)698:       IF (ALLOCATED(TWIST_THETA0)) DEALLOCATE(TWIST_THETA0)
697:       IF (ALLOCATED(TWIST_ATOMS)) DEALLOCATE(TWIST_ATOMS)699:       IF (ALLOCATED(TWIST_ATOMS)) DEALLOCATE(TWIST_ATOMS)
698: 700: 
699: !701: !
700: ! Should also deallocate any of the arrays below, which may have been702: ! Should also deallocate any of the arrays below, which may have been
701: ! allocated in keywords. This is just to make nagfmcheck happy, though! DJW703: ! allocated in keywords. This is just to make nagfmcheck happy, though! DJW
702: !704: !
703: !          ALLOCATE(DATOMGROUPAXIS(NDGROUPS,2))705: !          ALLOCATE(DATOMGROUPAXIS(NDGROUPS,2))


r31263/Makefile 2016-10-12 11:30:23.296774942 +0100 r31262/Makefile 2016-10-12 11:30:29.148853495 +0100
 51:         rotamer_move.o 51:         rotamer_move.o
 52: OBJS2 =        centre.o genrigid.o rotations.o finalio.o modconsts_trans_97.o modconsts.o dist.o tryexchange.o \ 52: OBJS2 =        centre.o genrigid.o rotations.o finalio.o modconsts_trans_97.o modconsts.o dist.o tryexchange.o \
 53:         io1.o keywords.o  main.o mc.o mcruns.o morse.o \ 53:         io1.o keywords.o  main.o mc.o mcruns.o morse.o \
 54:         potential.o quench.o dprand.o saveit.o seed.o \ 54:         potential.o quench.o dprand.o saveit.o seed.o \
 55:         sort.o sort2.o sort3.o sort4.o takestep.o mycpu_time.o trans.o \ 55:         sort.o sort2.o sort3.o sort4.o takestep.o mycpu_time.o trans.o \
 56:         finalq.o symmetry.o symmetrycsm.o ptgrp.o eigsrt.o SiSW.o taboo.o reseed.o newinertia.o supermc.o pgsym.o pgsym_mod.o\ 56:         finalq.o symmetry.o symmetrycsm.o ptgrp.o eigsrt.o SiSW.o taboo.o reseed.o newinertia.o supermc.o pgsym.o pgsym_mod.o\
 57:         tosifumi.o ortho.o compress.o mylbfgs.o mymylbfgs.o input.o ddfpmin.o dlnsrch.o cgmin.o linmin.o \ 57:         tosifumi.o ortho.o compress.o mylbfgs.o mymylbfgs.o input.o ddfpmin.o dlnsrch.o cgmin.o linmin.o \
 58:         brent.o mnbrak.o dbrent.o f1dim.o zwischen.o hsmove.o PachecoC60.o AT.o EAMLJ_sub.o \ 58:         brent.o mnbrak.o dbrent.o f1dim.o zwischen.o hsmove.o PachecoC60.o AT.o EAMLJ_sub.o \
 59:         Pbglue.o wenzel.o odesd.o capsid.o rigidfuncs.o tip.o pah.o strand.o \ 59:         Pbglue.o wenzel.o odesd.o capsid.o rigidfuncs.o tip.o pah.o strand.o \
 60:         SW.o qmod.o ljpbin.o fdm.o dftb.o ljpshift.o dzugutov.o ljcoulomb.o \ 60:         SW.o qmod.o ljpbin.o fdm.o dftb.o ljpshift.o dzugutov.o ljcoulomb.o \
 61:         enperms.o multiperm.o mie_field.o boxcentroid.o \ 61:         binary_id_swaps.o homoref.o homoref_addons.o enperms.o multiperm.o mie_field.o boxcentroid.o \
 62:         QALCSearch.o QALCS_surface.o QALCS_symmetry.o QALCS_symmetry_mod.o atom_label_swaps.o atom_label_flips.o stress.o \ 62:         QALCSearch.o QALCS_surface.o QALCS_symmetry.o QALCS_symmetry_mod.o atom_label_swaps.o atom_label_flips.o stress.o \
 63:         fd.o fedor.o welch.o glj_yukawa.o BGupta.o BLJcluster.o BLJcluster_nocut.o stock.o Farkas.o getorbits.o \ 63:         fd.o fedor.o welch.o glj_yukawa.o BGupta.o BLJcluster.o BLJcluster_nocut.o stock.o Farkas.o getorbits.o \
 64:         sc.o MSC.o MGupta.o MLJ.o Zetterling.o MSorig.o MSorigc.o MStrans.97.o convert.o \ 64:         sc.o MSC.o MGupta.o MLJ.o Zetterling.o MSorig.o MSorigc.o MStrans.97.o convert.o \
 65:         frausi.o p46merdiff.o g46merdiff.o lj.o modperm.o modf1com.o mododesd.o EAMal.o Alglue.o Mgglue.o \ 65:         frausi.o p46merdiff.o g46merdiff.o lj.o modperm.o modf1com.o mododesd.o EAMal.o Alglue.o Mgglue.o \
 66:         Gupta.o orient.o Natb.o sticky.o enumerate.o minperm.o minpermdist.o LB2.o \ 66:         Gupta.o orient.o Natb.o sticky.o enumerate.o minperm.o minpermdist.o LB2.o \
 67:         dgetrf.o dgetri.o reorient.o thomson.o Q4.o basinsampling.o tether.o tetherfunc.o BLN.o \ 67:         dgetrf.o dgetri.o reorient.o thomson.o Q4.o basinsampling.o tether.o tetherfunc.o BLN.o \
 68:         newmindist.o centrecom.o qorderparam_blj.o qorderparam_lj.o bspt.o GMINdump.o quad.o \ 68:         newmindist.o centrecom.o qorderparam_blj.o qorderparam_lj.o bspt.o GMINdump.o quad.o \
 69:         rigidbaa.o checkd.o capbin.o dbpg.o dbptd.o dmblmorse.o dmblpy.o gbcalamitic.o gbdiscotic.o gem.o gbddp.o linrod.o \ 69:         rigidbaa.o checkd.o capbin.o dbpg.o dbptd.o dmblmorse.o dmblpy.o gbcalamitic.o gbdiscotic.o gem.o gbddp.o linrod.o \
 70:         lwotp.o        newcapsid.o newpah.o msgayberne.o mstbin.o multstock.o paha.o pahw99.o pap.o ptstst.o papbin.o papjanus.o silane.o \ 70:         lwotp.o        newcapsid.o newpah.o msgayberne.o mstbin.o multstock.o paha.o pahw99.o pap.o ptstst.o papbin.o papjanus.o silane.o \
 71:         multpaha.o newtip.o patchy.o asaoos.o stockaa.o tetrahedra.o waterpdc.o waterpkz.o takestepmsgb.o dipolarmorse.o \ 71:         multpaha.o newtip.o patchy.o asaoos.o stockaa.o tetrahedra.o waterpdc.o waterpkz.o takestepmsgb.o dipolarmorse.o \
507: libcharmm.a: commons.o modcharmm.o modmxatms.o507: libcharmm.a: commons.o modcharmm.o modmxatms.o
508: libamber.a: commons.o modamber9.o porfuncs.o grouprotation.o508: libamber.a: commons.o modamber9.o porfuncs.o grouprotation.o
509: #${OBJS2}: ${OBJS1} 509: #${OBJS2}: ${OBJS1} 
510: 510: 
511: 511: 
512: Alglue.o:      commons.o512: Alglue.o:      commons.o
513: atomlists.o:   commons.o513: atomlists.o:   commons.o
514: atom_label_swaps.o: commons.o514: atom_label_swaps.o: commons.o
515: atom_label_flips.o: commons.o515: atom_label_flips.o: commons.o
516: boxcentroid.o: commons.o        516: boxcentroid.o: commons.o        
 517: binary_id_swaps.o: commons.o
517: BLJcluster.o:  commons.o518: BLJcluster.o:  commons.o
518: GLJ.o:  commons.o519: GLJ.o:  commons.o
519: BLJcluster_nocut.o: commons.o520: BLJcluster_nocut.o: commons.o
520: BGupta.o:      commons.o521: BGupta.o:      commons.o
521: detsym.o:      commons.o522: detsym.o:      commons.o
 523: homoref.o:     commons.o
 524: homoref_addons.o: commons.o
522: enperms.o:     commons.o525: enperms.o:     commons.o
523: multiperm.o:   commons.o porfuncs.o526: multiperm.o:   commons.o porfuncs.o
524: EAMLJ_sub.o:   commons.o527: EAMLJ_sub.o:   commons.o
525: EAMal.o:       commons.o528: EAMal.o:       commons.o
526: Gupta.o:       commons.o529: Gupta.o:       commons.o
527: glj_yukawa.o:  commons.o530: glj_yukawa.o:  commons.o
528: MSorig.o:      commons.o modconsts.o dist.o531: MSorig.o:      commons.o modconsts.o dist.o
529: MSorigc.o:     commons.o modconsts.o dist.o532: MSorigc.o:     commons.o modconsts.o dist.o
530: MStrans.97.o:  commons.o modconsts_trans_97.o dist.o533: MStrans.97.o:  commons.o modconsts_trans_97.o dist.o
531: Mgglue.o:      commons.o534: Mgglue.o:      commons.o


r31263/mc.F 2016-10-12 11:30:26.660820096 +0100 r31262/mc.F 2016-10-12 11:30:32.584899614 +0100
994: ! MARK'S MIXED CLUSTER STEP TAKING994: ! MARK'S MIXED CLUSTER STEP TAKING
995: !995: !
996: !     MAM1000 >  Charged/neutral particle swaps for LJ+Coulomb mixed clusters996: !     MAM1000 >  Charged/neutral particle swaps for LJ+Coulomb mixed clusters
997:                ELSE IF (LJCOULT) THEN997:                ELSE IF (LJCOULT) THEN
998:                   RANDOM=DPRAND()998:                   RANDOM=DPRAND()
999:                   IF (RANDOM < COULSWAP) THEN999:                   IF (RANDOM < COULSWAP) THEN
1000:                      CALL TAKESTEPLJC(JP)1000:                      CALL TAKESTEPLJC(JP)
1001:                      MCTEMP = COULTEMP1001:                      MCTEMP = COULTEMP
1002:                   ELSE1002:                   ELSE
1003:                      CALL TAKESTEP(JP)1003:                      CALL TAKESTEP(JP)
1004:                   END IF                  1004:                   END IF
 1005:                   
1005: !1006: !
 1007: ! ds656> Identity swaps for binary systems. If BASWAP is FALSE or
 1008: !        J1 <= NWAIT then the usual TAKESTEP will be called further down.
 1009: 
 1010:                ELSE IF (BASWAP .AND. (J1 > BASWAP_NWAIT)) THEN
 1011:                   !
 1012:                   RANDOM=DPRAND()
 1013:                   !IF(BASWAP_FRAC > 0.0D0) THEN ! no longer needed!!!
 1014:                   !   RANDOM=RANDOM/DEXP( (QMIN(1)-EPREV(JP))/BASWAP_TEMP )
 1015:                   !ENDIF 
 1016:                   !
 1017:                   IF(BASWAP_NMAX > 0 .AND. RANDOM < BASWAP_FRAC ) THEN
 1018:                      !
 1019:                      NSWAPS = 1 + INT(BASWAP_NMAX*DPRAND())
 1020:                      !     
 1021:                      DO J4=1,NSWAPS
 1022:                         CALL ABSWAP(JP, DE1, DE2)
 1023:                      ENDDO
 1024:                      !
 1025:                      IF(BASWAPTEST) THEN
 1026:                         !
 1027:                         IF(BGUPTAT) THEN
 1028:                            CALL BGUPTA(COORDS,GRAD,EASWAP,.FALSE.)
 1029:                         ELSEIF(BLJCLUSTER_NOCUT) THEN
 1030:                            CALL BLJ_CLUST(COORDS,GRAD,EASWAP,.FALSE.)
 1031:                         ENDIF
 1032:                         !
 1033:                      ENDIF
 1034:                      !
 1035:                   ELSEIF(BASWAP_NMAX==0 .AND. RANDOM < BASWAP_FRAC) THEN
 1036:                      !
 1037:                      NSWAPS = 1
 1038:                      !
 1039:                      !CALL BGUPTA(COORDS,GRAD,EBSWAP,.FALSE.)
 1040:                      CALL ABSWAP2(JP, DE1, DE2)
 1041:                      !
 1042:                      IF(BASWAPTEST .AND. BGUPTAT) THEN
 1043:                         CALL BGUPTA(COORDS,GRAD,EASWAP,.FALSE.)
 1044:                      ELSEIF(BASWAPTEST .AND. BLJCLUSTER_NOCUT) THEN
 1045:                         CALL BLJ_CLUST(COORDS,GRAD,EASWAP,.FALSE.)
 1046:                      ENDIF
 1047:                      !     
 1048:                   ELSE
 1049:                      CALL TAKESTEP(JP)
 1050:                   END IF
 1051: ! <ds656
1006: !               ELSE IF (LWOTPT) THEN1052: !               ELSE IF (LWOTPT) THEN
1007: !                  CALL TAKESTEPLWOTP(JP)1053: !                  CALL TAKESTEPLWOTP(JP)
1008: 1054: 
1009:                ELSE IF (GBT .OR. GBDT .OR. GBDPT) THEN1055:                ELSE IF (GBT .OR. GBDT .OR. GBDPT) THEN
1010:                   CALL TKSTDCELPSD(JP)1056:                   CALL TKSTDCELPSD(JP)
1011:       1057:       
1012:                ELSE IF (MSGBT) THEN1058:                ELSE IF (MSGBT) THEN
1013:                   CALL TAKESTEPMSGB(JP) 1059:                   CALL TAKESTEPMSGB(JP) 
1014: 1060: 
1015:                ELSE IF (SANDBOXT) THEN1061:                ELSE IF (SANDBOXT) THEN
1297: !  These coordinates are overwritten if we try to call MPI_IPROBE1343: !  These coordinates are overwritten if we try to call MPI_IPROBE
1298: !1344: !
1299: !                 WRITE(MYUNIT,'(I6)') NATOMS1345: !                 WRITE(MYUNIT,'(I6)') NATOMS
1300: !                 WRITE(MYUNIT,'(A,G20.10)') 'eprev=',EPREV1346: !                 WRITE(MYUNIT,'(A,G20.10)') 'eprev=',EPREV
1301: !                 WRITE(MYUNIT,'(A,3G20.10)') ('LA ',COORDS(3*(J3-1)+1:3*(J3-1)+3,JP),J3=1,NATOMS)1347: !                 WRITE(MYUNIT,'(A,3G20.10)') ('LA ',COORDS(3*(J3-1)+1:3*(J3-1)+3,JP),J3=1,NATOMS)
1302:                   IF (USERPOTT.AND.GROUPROTT.AND.DOGROUPROT) CALL GROUPROTSTEP(JP)1348:                   IF (USERPOTT.AND.GROUPROTT.AND.DOGROUPROT) CALL GROUPROTSTEP(JP)
1303: !ab2111> dihedralrotation with user-defined potential1349: !ab2111> dihedralrotation with user-defined potential
1304:                   IF (USERPOTT.AND.DIHEDRALROTT.AND.DODIHEDRALROT) CALL DIHEDRALROTSTEP(DUMMY1,JP,DUMMY2,DUMMY3)1350:                   IF (USERPOTT.AND.DIHEDRALROTT.AND.DODIHEDRALROT) CALL DIHEDRALROTSTEP(DUMMY1,JP,DUMMY2,DUMMY3)
1305:                   CALL TAKESTEP(JP)1351:                   CALL TAKESTEP(JP)
1306: !ds656>1352: !ds656>
 1353:                   IF (RANDPERMT .AND. MOD(J1,RANDPERM_STEP)==0) CALL RANDPERM(JP)
1307:                   IF (RANDMULTIPERMT.AND.MOD(J1,RANDMULTIPERM_STEP)==0) CALL RANDMULTIPERM(JP,1)1354:                   IF (RANDMULTIPERMT.AND.MOD(J1,RANDMULTIPERM_STEP)==0) CALL RANDMULTIPERM(JP,1)
1308: 1355: 
1309:                ENDIF1356:                ENDIF
1310: ! jdf43> 1357: ! jdf43> 
1311:                IF (DDMT) THEN1358:                IF (DDMT) THEN
1312:                   DO I2=1,3*NATOMS1359:                   DO I2=1,3*NATOMS
1313:                      CALL DDMCONDENSE(COORDS)1360:                      CALL DDMCONDENSE(COORDS)
1314:                   ENDDO1361:                   ENDDO
1315:                ENDIF1362:                ENDIF
1316: ! Restore atom coordinates if atom is FROZEN or DONTMOVE as long as1363: ! Restore atom coordinates if atom is FROZEN or DONTMOVE as long as
1426: !       Markov minimum, but not too far...1473: !       Markov minimum, but not too far...
1427:                IF(SAWTOOTH) THEN1474:                IF(SAWTOOTH) THEN
1428:                   ! Energy criterion alone may not be adequate...1475:                   ! Energy criterion alone may not be adequate...
1429:                   IF(ABS(POTEL-EPREV(JP)) < ECONV) THEN1476:                   IF(ABS(POTEL-EPREV(JP)) < ECONV) THEN
1430:                      STEP(JP)=SAWTOOTH_SFAC*STEP(JP)1477:                      STEP(JP)=SAWTOOTH_SFAC*STEP(JP)
1431:                      GOTO 23 ! Repeat random step1478:                      GOTO 23 ! Repeat random step
1432:                   ELSE1479:                   ELSE
1433:                      STEP(JP)=SAWTOOTH_SFAC2*STEP(JP)1480:                      STEP(JP)=SAWTOOTH_SFAC2*STEP(JP)
1434:                   ENDIF1481:                   ENDIF
1435:                   WRITE(MYUNIT,'(A,F10.5)') 'mc> STEP=', STEP(JP)1482:                   WRITE(MYUNIT,'(A,F10.5)') 'mc> STEP=', STEP(JP)
 1483:                ENDIF
 1484: !ds656> print some energies to test BASWAP
 1485:                IF(BASWAPTEST) THEN
 1486:                   !CALL BGUPTA(COORDS,GRAD,EASWAPQ,.FALSE.)
 1487:                   IF(NSWAPS == 1) THEN
 1488:                      WRITE(MYUNIT, '(2(A,1X,E12.6,1X))') 
 1489:      1                    "mc> 1 swap: de_A=",DE1,"and de_B=",DE2
 1490:                   ENDIF
 1491:                   WRITE(MYUNIT, '(2(A,1X,E12.6,1X),A)') 
 1492:      2                 "mc> BA swaps: dE=", EASWAP-EPREV(JP), "before and dE=",
 1493:      3                 POTEL-EPREV(JP), "after quench."
 1494:                   !WRITE(MYUNIT, '(A,3(1X,F15.8))') "mc> ds656:", EPREV(JP)-EBSWAP, EASWAP, EASWAPQ-POTEL 
 1495:                   NSWAPS=0
1436:                ENDIF               1496:                ENDIF               
 1497: !<ds656        
 1498:                
1437: !1499: !
1438: !  Output1500: !  Output
1439: !1501: !
1440: !fh301>{{{1502: !fh301>{{{
1441:                IF (CHEMSHIFT2) THEN1503:                IF (CHEMSHIFT2) THEN
1442:                   CHEMSHIFT=.TRUE.1504:                   CHEMSHIFT=.TRUE.
1443: !     CHEMSHIFTITER=01505: !     CHEMSHIFTITER=0
1444:                ENDIF1506:                ENDIF
1445:                IF (MOD(J1-1,PRTFRQ).EQ.0) THEN !<start of looong block of print statements1507:                IF (MOD(J1-1,PRTFRQ).EQ.0) THEN !<start of looong block of print statements
1446:                   IF (CHEMSHIFT) THEN1508:                   IF (CHEMSHIFT) THEN
1548:      1              (MOD(J1,LSWAPS_N1)==0) ) THEN1610:      1              (MOD(J1,LSWAPS_N1)==0) ) THEN
1549:                   CALL USWAPS(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)1611:                   CALL USWAPS(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
1550:                ENDIF1612:                ENDIF
1551: 1613: 
1552:                IF ( LFLIPST .AND.1614:                IF ( LFLIPST .AND.
1553:      1              (MOD(J1,LFLIPS_N1)==0) ) THEN1615:      1              (MOD(J1,LFLIPS_N1)==0) ) THEN
1554:                   CALL FLIPSEQ(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)1616:                   CALL FLIPSEQ(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
1555:                ENDIF1617:                ENDIF
1556: 1618: 
1557:                IF ((QDONE.EQ.1).AND.(.NOT.SYMMETRIZE)) THEN1619:                IF ((QDONE.EQ.1).AND.(.NOT.SYMMETRIZE)) THEN
 1620:                   IF (HOMOREF_BHT) CALL HOMOREF_BH(JP, ITERATIONS, TIME, BRUN, QDONE, SCREENC)
 1621:                   IF (HOMOREFT) CALL HOMOREF(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)                 
1558:                   IF (QALCST.OR.QALCS_SURFT.OR.QALCS_SYMT) CALL QALCS(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)1622:                   IF (QALCST.OR.QALCS_SURFT.OR.QALCS_SYMT) CALL QALCS(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
1559:                ENDIF1623:                ENDIF
1560: !<ds6561624: !<ds656
1561:                1625:                
1562: !  RMS compared to reference structure 'compare'1626: !  RMS compared to reference structure 'compare'
1563:           IF (RMST.AND.CHRMMT) THEN1627:           IF (RMST.AND.CHRMMT) THEN
1564:              CALL CHRMS(JP,RMSD)1628:              CALL CHRMS(JP,RMSD)
1565:              IF (DEBUG) WRITE(MYUNIT,'(A,F15.5)')'RMSD = ',RMSD1629:              IF (DEBUG) WRITE(MYUNIT,'(A,F15.5)')'RMSD = ',RMSD
1566:              IF (RMSD.LE.RMSLIMIT) CALL SAVERMS(JP,POTEL,RMSD)1630:              IF (RMSD.LE.RMSLIMIT) CALL SAVERMS(JP,POTEL,RMSD)
1567: !                NRMS=NRMS+11631: !                NRMS=NRMS+1


r31263/mc_gbh.F90 2016-10-12 11:30:26.916823535 +0100 r31262/mc_gbh.F90 2016-10-12 11:30:32.840903048 +0100
  1: !   GMIN: A program for finding global minima  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/mc_gbh.F90' in revision 31262
  2: !   Copyright (C) 1999-2006 David J. Wales 
  3: !   This file is part of GMIN and was writtedn by ds656.    
  4: ! 
  5: !   GMIN is free software; you can redistribute it and/or modify 
  6: !   it under the terms of the GNU General Public License as published by 
  7: !   the Free Software Foundation; either version 2 of the License, or    
  8: !   (at your option) any later version. 
  9: ! 
 10: !   GMIN is distributed in the hope that it will be useful,    
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: !   GNU General Public License for more details. 
 14: !                                                
 15: !   You should have received a copy of the GNU General Public License 
 16: !   along with this program; if not, write to the Free Software     
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 18: ! 
 19: SUBROUTINE MC_GBH(NSTEPS) 
 20:   ! 
 21:   USE PORFUNCS 
 22:   USE COMMONS, ONLY : MYNODE, MYUNIT, DEBUG, NATOMSALLOC, NQ, & 
 23:        TSTART, COORDS, LABELS, ECONV, BOXCENTROIDT,RANDMULTIPERMT, & 
 24:        NPAR_GBH, TEMP, RMS, TARGET, HIT, QALCST, QALCS_NBRHD, & 
 25:        QALCSV 
 26:   ! 
 27:   IMPLICIT NONE 
 28:   ! 
 29: #ifdef MPI  
 30:   INCLUDE 'mpif.h' 
 31:   INTEGER MPIERR 
 32: #endif 
 33:   ! 
 34:   INTEGER, INTENT(IN) :: NSTEPS 
 35:   ! 
 36:   LOGICAL :: ACCEPT 
 37:   INTEGER :: I, ITRIAL, ITERNS, BRUN, QDONE, NQTOT, IPROC, IPROCLO 
 38:   DOUBLE PRECISION :: TIME, SCREENC(3*NATOMSALLOC), POTEL, & 
 39:        POTEL_LIST(NPAR_GBH), POTELO, ELOWEST, R, DPRAND 
 40:   ! 
 41:   COMMON /MYPOT/ POTEL 
 42:   ! 
 43:   WRITE(MYUNIT, '(A)')  'mc_gbh> Calculating initial energy' 
 44:   !WRITE(MYUNIT, *)  'mc_gbh> NATOMSALLOC=', NATOMSALLOC 
 45:   CALL FLUSH(MYUNIT) 
 46:   CALL QUENCH(.FALSE.,1,ITERNS,TIME,BRUN,QDONE,SCREENC) 
 47:   NQ(1) = 0 
 48:   WRITE(MYUNIT,& 
 49:        '(A,I10,A,G20.10,A,I5,A,G12.5,A,F10.1)') & 
 50:        'Qu ',NQ(1),' E= ',POTEL,' steps= ',ITERNS,' RMS= ',RMS,& 
 51:        ' t= ',TIME-TSTART 
 52:   CALL FLUSH(MYUNIT) 
 53:   ! 
 54:   WRITE(MYUNIT, '(A,I6)') & 
 55:        'mc_gbh> Starting GBH loop of length ',NSTEPS 
 56:   ! 
 57:   gbh_loop: DO ITRIAL=1,NSTEPS 
 58:      ! 
 59:      IF(QALCSV) WRITE(MYUNIT, '(A,I10)') & 
 60:           'mc_gbh> Starting iteration ', ITRIAL      
 61:      POTELO=POTEL 
 62:      ! 
 63:      ! --- Stochastic cartesian move ------- 
 64:      CALL TAKESTEP(1) 
 65:      !IF (RANDMULTIPERMT.AND.MOD(J1,RANDMULTIPERM_STEP)==0) & 
 66:      !     CALL RANDMULTIPERM(1) ! re-write this routine!!!!! 
 67:      !IF(BOXCENTROIDT) CALL BOXCENTROID(COORDS(:,1)) 
 68:      ! 
 69:      ! Quench perturbed state 
 70:      NQ(1) = NQ(1) + 1 
 71:      CALL QUENCH(.FALSE.,1,ITERNS,TIME,BRUN,QDONE,SCREENC) 
 72:      WRITE(MYUNIT,& 
 73:           '(A,I10,A,G20.10,A,I5,A,G12.5,A,F10.1)') & 
 74:           'Qu ',NQ(1),' E= ',POTEL,' steps= ',ITERNS,' RMS= ',& 
 75:           RMS,' t= ',TIME-TSTART 
 76:      CALL FLUSH(MYUNIT) 
 77:      ! 
 78: #ifdef MPI 
 79:      ! Gather all parallel quench energies on master node. 
 80:      CALL MPI_GATHER(POTEL, 1, MPI_DOUBLE_PRECISION, & 
 81:           POTEL_LIST, 1, MPI_DOUBLE_PRECISION, 0, & 
 82:           MPI_COMM_WORLD, MPIERR) 
 83: #else 
 84:      POTEL_LIST(1) = POTEL 
 85: #endif 
 86:      ! 
 87:      IF(MYNODE==0) THEN 
 88:         ! Find lowest energy in the gathered list. 
 89:         ELOWEST=MINVAL(POTEL_LIST) 
 90:         IPROCLO=MINLOC(POTEL_LIST, 1)-1 
 91:         ! Choose a IPROC with Boltzmann probability. 
 92:         CALL CHOOSE_FROM_LIST(NPAR_GBH,POTEL_LIST,IPROC) 
 93:         IPROC=IPROC-1 
 94:      ENDIF 
 95:      ! 
 96: #ifdef MPI  
 97:      ! Broadcast IPROC from master 
 98:      CALL MPI_BCAST(IPROC,1,MPI_INTEGER,& 
 99:           0,MPI_COMM_WORLD,MPIERR) 
100:      ! Broadcast state from IPROC 
101:      CALL MPI_BCAST(COORDS(:,1),3*NATOMSALLOC,& 
102:           MPI_DOUBLE_PRECISION,IPROC,MPI_COMM_WORLD,MPIERR) 
103:      CALL MPI_BCAST(LABELS(:,1),NATOMSALLOC,MPI_INTEGER,& 
104:           IPROC,MPI_COMM_WORLD,MPIERR) 
105:      CALL MPI_BCAST(POTEL,1,MPI_DOUBLE_PRECISION,& 
106:           IPROC,MPI_COMM_WORLD,MPIERR) 
107:      WRITE(MYUNIT, '(A,F20.10,A,I6)') & 
108:           'mc_gbh> Use E= ', POTEL,' from node ',IPROC+1 
109:      ! Broadcast HIT from IPROCLO if targetting... 
110:      IF(TARGET) THEN 
111:         ! Broadcast IPROCLO from master 
112:         CALL MPI_BCAST(IPROCLO,1,MPI_INTEGER,& 
113:              0,MPI_COMM_WORLD,MPIERR) 
114:         CALL MPI_BCAST(HIT,1,MPI_LOGICAL,& 
115:              IPROCLO,MPI_COMM_WORLD,MPIERR)         
116:      ENDIF 
117: #endif                      
118:      ! 
119:      IF(HIT) THEN 
120:         WRITE(MYUNIT,'(A,I3,A,I6)') & 
121:              'mc_gbh> Target hit stochastically by node ',IPROCLO+1,& 
122:              ' on trial ',ITRIAL 
123:         EXIT gbh_loop 
124:      ENDIF 
125:      ! 
126:      ! --- Deterministic refinement -------------------------- 
127:      IF(.NOT.HIT .AND. QALCST .AND. QALCS_NBRHD>0 & 
128:           .AND. ABS(POTEL-POTELO)>ECONV) THEN 
129:         ! perform a variable neighbourhood search (in parallel) 
130:         CALL QALCS_PARALLEL(ITERNS, TIME, BRUN, QDONE, SCREENC) 
131:      ENDIF 
132:      ! 
133:      IF(HIT) EXIT gbh_loop 
134:      ! 
135:      ! 
136:   ENDDO gbh_loop 
137:   ! 
138:   RETURN 
139:   ! 
140: END SUBROUTINE MC_GBH 
141:  
142: SUBROUTINE CHOOSE_FROM_LIST(N,VALUES,I) 
143:   ! 
144:   ! Choose an element of VALUES with Boltzmann probability 
145:   ! 
146:   USE COMMONS, ONLY : TEMP 
147:   ! 
148:   IMPLICIT NONE 
149:   ! 
150:   INTEGER, INTENT(IN) :: N 
151:   DOUBLE PRECISION, INTENT(IN) :: VALUES(N) 
152:   INTEGER, INTENT(OUT) :: I 
153:   ! 
154:   INTEGER :: J 
155:   DOUBLE PRECISION :: PSUM(N), X, DPRAND 
156:   ! 
157:   X=0.0D0 ! initialise total sum 
158:   DO J=1,N  
159:      X = X + DEXP(-VALUES(J)/TEMP(1)) 
160:      PSUM(J) = X ! store partial sum 
161:   ENDDO 
162:   ! 
163:   X=X*DPRAND() 
164:   ! 
165:   I=1 
166:   DO WHILE (X > PSUM(I)) 
167:      I=I+1 
168:   ENDDO 
169:   ! 
170:   RETURN 
171:   ! 
172: END SUBROUTINE CHOOSE_FROM_LIST 


r31263/mcruns.f90 2016-10-12 11:30:27.168826915 +0100 r31262/mcruns.f90 2016-10-12 11:30:33.156907294 +0100
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13: !   GNU General Public License for more details. 13: !   GNU General Public License for more details.
 14: ! 14: !
 15: !   You should have received a copy of the GNU General Public License 15: !   You should have received a copy of the GNU General Public License
 16: !   along with this program; if not, write to the Free Software 16: !   along with this program; if not, write to the Free Software
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19: SUBROUTINE MCRUNS(SCREENC) 19: SUBROUTINE MCRUNS(SCREENC)
 20:    USE COMMONS, ONLY: NATOMS, PTMC, BSPT, ONE_ATOM_TAKESTEP, CHECKDT, COORDS, & 20:    USE COMMONS, ONLY: NATOMS, PTMC, BSPT, ONE_ATOM_TAKESTEP, CHECKDT, COORDS, &
 21:                       NRUNS, MCSTEPS, TFAC, RELAXFQ, NSAVE, GBHT 21:                       NRUNS, MCSTEPS, TFAC, RELAXFQ, NSAVE
 22:    USE GENRIGID, ONLY: RIGIDINIT 22:    USE GENRIGID, ONLY: RIGIDINIT
 23:    USE QMODULE, ONLY: QMINP 23:    USE QMODULE, ONLY: QMINP
 24:    USE PREC, ONLY: INT32, REAL64 24:    USE PREC, ONLY: INT32, REAL64
 25:    IMPLICIT NONE 25:    IMPLICIT NONE
 26: ! Arguments (what intent for SCREENC?) 26: ! Arguments (what intent for SCREENC?)
 27:    REAL(REAL64)               :: SCREENC(3*NATOMS) 27:    REAL(REAL64)               :: SCREENC(3*NATOMS)
 28: ! Variables 28: ! Variables
 29:    LOGICAL                    :: LOPEN 29:    LOGICAL                    :: LOPEN
 30:    REAL(REAL64), ALLOCATABLE  :: QMINPTMP(:, :) 30:    REAL(REAL64), ALLOCATABLE  :: QMINPTMP(:, :)
 31:    INTEGER(INT32)             :: I 31:    INTEGER(INT32)             :: I
 32:  32: 
 33:    IF (PTMC .OR. BSPT) THEN 33:    IF (PTMC .OR. BSPT) THEN
 34:       IF (ONE_ATOM_TAKESTEP) THEN 34:       IF (ONE_ATOM_TAKESTEP) THEN
 35:          CALL PTMC_ONE_ATOM() 35:          CALL PTMC_ONE_ATOM()
 36:       ELSE 36:       ELSE
 37:          CALL PTBASINSAMPLING() 37:          CALL PTBASINSAMPLING()
 38:       END IF 38:       END IF
 39:       RETURN 39:       RETURN
 40:    END IF 40:    END IF
 41:  41: 
 42:    IF (GBHT) THEN ! ds656> Generalised basin-hopping run. 
 43:       CALL MC_GBH(MCSTEPS(1)) 
 44:       CALL FINALQ() 
 45:       CALL FINALIO() 
 46:       RETURN 
 47:    ENDIF 
 48:  
 49: !   IF (MFETT) THEN 42: !   IF (MFETT) THEN
 50: !      CALL MFET() 43: !      CALL MFET()
 51: !      CALL FINALIO() 44: !      CALL FINALIO()
 52: !      RETURN 45: !      RETURN
 53: !   ENDIF 46: !   ENDIF
 54:  47: 
 55:    INQUIRE(UNIT=1, OPENED=LOPEN) 48:    INQUIRE(UNIT=1, OPENED=LOPEN)
 56:    IF (LOPEN) THEN 49:    IF (LOPEN) THEN
 57:       WRITE(*,'(A,I2,A)') 'mcruns> ERROR *** Unit ', 1, ' is not free ' 50:       WRITE(*,'(A,I2,A)') 'mcruns> ERROR *** Unit ', 1, ' is not free '
 58:       STOP 51:       STOP


r31263/MLJ.f90 2016-10-12 11:30:23.024771289 +0100 r31262/MLJ.f90 2016-10-12 11:30:28.884849907 +0100
 19: !  =============================================================== 19: !  ===============================================================
 20: !  Multicomponent Lennard-Jones without a cutoff. 20: !  Multicomponent Lennard-Jones without a cutoff.
 21: !  Assumed reduced units with SIGMA_AA = EPSILON_AA = 1. 21: !  Assumed reduced units with SIGMA_AA = EPSILON_AA = 1.
 22: !  Per-atom energy is stored in array VT(NATOMS). 22: !  Per-atom energy is stored in array VT(NATOMS).
 23: ! 23: !
 24: !  ds656> 6/12/2014 24: !  ds656> 6/12/2014
 25: !  =============================================================== 25: !  ===============================================================
 26: ! 26: !
 27: SUBROUTINE MLJ(X,GRAD,POT,GRADT,HESST,STRESST) 27: SUBROUTINE MLJ(X,GRAD,POT,GRADT,HESST,STRESST)
 28:   ! 28:   !
 29:   USE COMMONS, ONLY : NATOMS,VT,LABELS,NSPECIES,STRESS,MYNODE,GBHT 29:   USE COMMONS, ONLY : NATOMS,VT,LABELS,NSPECIES,STRESS,MYNODE
 30:   USE MODHESS, ONLY : HESS 30:   USE MODHESS, ONLY : HESS
 31:   USE POT_PARAMS, ONLY : MLJ_EPS, MLJ_SIG 31:   USE POT_PARAMS, ONLY : MLJ_EPS, MLJ_SIG
 32:   ! 32:   !
 33:   IMPLICIT NONE 33:   IMPLICIT NONE
 34:   ! 34:   !
 35:   LOGICAL, INTENT (IN) :: GRADT,HESST,STRESST 35:   LOGICAL, INTENT (IN) :: GRADT,HESST,STRESST
 36:   DOUBLE PRECISION, INTENT (IN) :: X(3*NATOMS) 36:   DOUBLE PRECISION, INTENT (IN) :: X(3*NATOMS)
 37:   DOUBLE PRECISION, INTENT (INOUT) :: POT, GRAD(3*NATOMS) 37:   DOUBLE PRECISION, INTENT (INOUT) :: POT, GRAD(3*NATOMS)
 38:   ! 38:   !
 39:   INTEGER :: I,J,J1,J13,J2,J23,T1,T2,IP 39:   INTEGER :: I,J,J1,J13,J2,J23,T1,T2,IP
 47:   POT = 0.0D0 47:   POT = 0.0D0
 48:   IF(GRADT) GRAD(:) = 0.0D0 48:   IF(GRADT) GRAD(:) = 0.0D0
 49:   IF(HESST) HESS(:,:) = 0.0D0 49:   IF(HESST) HESS(:,:) = 0.0D0
 50:   IF(STRESST) STRESS(:,:,:) = 0.0D0 50:   IF(STRESST) STRESS(:,:,:) = 0.0D0
 51:   ! 51:   !
 52:   ! This can be precomputed elsewhere... 52:   ! This can be precomputed elsewhere...
 53:   SIG2(:,:) = MLJ_SIG(:,:)*MLJ_SIG(:,:) 53:   SIG2(:,:) = MLJ_SIG(:,:)*MLJ_SIG(:,:)
 54:   EPS4(:,:) = 4.0D0*MLJ_EPS(:,:) 54:   EPS4(:,:) = 4.0D0*MLJ_EPS(:,:)
 55:   EPS24(:,:) = 24.0D0*MLJ_EPS(:,:) 55:   EPS24(:,:) = 24.0D0*MLJ_EPS(:,:)
 56:   ! 56:   !
 57:   IF(GBHT) THEN 57:   IP=MYNODE+1
 58:      IP=1 
 59:   ELSE 
 60:      IP=MYNODE+1 
 61:   ENDIF 
 62:   ! 58:   !
 63:   ! Double loop over atom types 59:   ! Double loop over atom types
 64:   DO J1=1,NATOMS-1 60:   DO J1=1,NATOMS-1
 65:      J13=3*(J1-1) 61:      J13=3*(J1-1)
 66:      T1=LABELS(J1,IP) 62:      T1=LABELS(J1,IP)
 67:      DO J2=J1+1,NATOMS 63:      DO J2=J1+1,NATOMS
 68:         J23 = 3*(J2-1) 64:         J23 = 3*(J2-1)
 69:         T2=LABELS(J2,IP) 65:         T2=LABELS(J2,IP)
 70:         ! 66:         !
 71:         ! 67:         !


r31263/multiperm.f90 2016-10-12 11:30:27.428830406 +0100 r31262/multiperm.f90 2016-10-12 11:30:33.412910741 +0100
101:      !101:      !
102:      COORDS(:,JP) = X0(:) ! re-initialise the coordinates102:      COORDS(:,JP) = X0(:) ! re-initialise the coordinates
103:      CALL MULTIPERM_NEXT(NATOMS,LIST,MORE) ! permutation step103:      CALL MULTIPERM_NEXT(NATOMS,LIST,MORE) ! permutation step
104:      CALL SET_ATOMLISTS(LIST,1) ! update atom lists104:      CALL SET_ATOMLISTS(LIST,1) ! update atom lists
105:      !105:      !
106:   ENDDO106:   ENDDO
107:   ! ============================================================107:   ! ============================================================
108:   !108:   !
109:   WRITE(MYUNIT, '(A,I10)')  &109:   WRITE(MYUNIT, '(A,I10)')  &
110:        'multiperm> Finished with permutation count=', K110:        'multiperm> Finished with permutation count=', K
111:   CALL FINALQ()  111:   CALL FINALQ  
112:   CALL FINALIO()112:   CALL FINALIO 
113:   !113:   !
114:   RETURN114:   RETURN
115:   !115:   !
116: END SUBROUTINE MULTIPERM116: END SUBROUTINE MULTIPERM
117: !117: !
118: ! ================================================================118: ! ================================================================
119: ! ================================================================119: ! ================================================================
120: ! ds656> All the routines that follow have been taken from120: ! ds656> All the routines that follow have been taken from
121: ! http://people.sc.fsu.edu/~jburkardt/f_src/subset/subset.f90121: ! http://people.sc.fsu.edu/~jburkardt/f_src/subset/subset.f90
122: ! ================================================================122: ! ================================================================


r31263/QALCS_parallel.F90 2016-10-12 11:30:23.564778534 +0100 r31262/QALCS_parallel.F90 2016-10-12 11:30:29.408856990 +0100
  1: SUBROUTINE QALCS_PARALLEL(ITER,TIME,BRUN,QDONE,SCREENC)  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/QALCS_parallel.F90' in revision 31262
  2:   ! 
  3:   USE COMMONS, ONLY : MYUNIT, NATOMS 
  4:   ! 
  5:   IMPLICIT NONE 
  6:   ! 
  7:   ! Passed variables for QUENCH 
  8:   INTEGER, INTENT(INOUT) :: ITER, BRUN, QDONE 
  9:   DOUBLE PRECISION, INTENT(INOUT) :: TIME, SCREENC(3*NATOMS) 
 10:   ! 
 11:   LOGICAL :: SWAP 
 12:   INTEGER :: I 
 13:   DOUBLE PRECISION :: POTEL 
 14:   COMMON /MYPOT/ POTEL   
 15:   ! 
 16:   I=0 
 17:   SWAP=.TRUE. 
 18:   DO WHILE(SWAP) 
 19:      CALL SCAN_ALL_SWAPS_PARALLEL(ITER,TIME,BRUN,QDONE,SCREENC,SWAP) 
 20:      IF(SWAP) WRITE(MYUNIT,'(A, G20.10)') & 
 21:           'QALCS_parallel> New E= ', POTEL 
 22:      I=I+1 
 23:   ENDDO 
 24:   WRITE(MYUNIT,'(A,I6,A)') 'QALCS_parallel> Converged in ', I, & 
 25:        ' swaps.' 
 26:   ! 
 27:   RETURN 
 28:   ! 
 29: END SUBROUTINE QALCS_PARALLEL 
 30:  
 31: !SUBROUTINE EXE_BEST_SWAP(ITER,TIME,BRUN,QDONE,SCREENC,SWAP) 
 32:   ! 
 33:   ! 
 34: !END SUBROUTINE EXE_BEST_SWAP 
 35:  
 36: SUBROUTINE SCAN_ALL_SWAPS_PARALLEL(ITER,TIME,BRUN,QDONE,SCREENC,SWAP) 
 37:   ! 
 38:   USE COMMONS, ONLY : NATOMS,COORDS,COORDSO,LABELS,LABELSO, & 
 39:        MYUNIT, MYNODE, QALCS_NBRHD, NPAR_GBH, ECONV, NQ, RMS, & 
 40:        QALCSV, TARGET, HIT 
 41:   ! 
 42:   IMPLICIT NONE 
 43:   ! 
 44: #ifdef MPI 
 45:   INCLUDE 'mpif.h' 
 46:   INTEGER MPIERR 
 47: #endif 
 48:   ! 
 49:   ! Passed variables for QUENCH 
 50:   INTEGER, INTENT(INOUT) :: ITER, BRUN, QDONE 
 51:   DOUBLE PRECISION, INTENT(INOUT) :: TIME, SCREENC(3*NATOMS) 
 52:   ! Indicator of whether an improvement has been made 
 53:   LOGICAL, INTENT(OUT) :: SWAP 
 54:   ! 
 55:   INTEGER :: I, J, K, NBRHOOD(2*QALCS_NBRHD), NPPN, MYI1, MYI2, & 
 56:        LMIN(NATOMS), IPROC 
 57:   DOUBLE PRECISION :: POTEL, POTELO, EMIN, XMIN(3*NATOMS), & 
 58:        POTEL_LIST(NPAR_GBH), ELOWEST 
 59:   ! 
 60:   ! Energy of COORDS from last quench. Common block in QUENCH. 
 61:   COMMON /MYPOT/ POTEL   
 62:   ! 
 63:   IF(QALCSV) WRITE(MYUNIT,'(A)') & 
 64:        'scan_all_swaps_parallel> Building neighbourhood list...' 
 65:   ! First build the list of atom pairs to be swapped 
 66:   K=1 
 67:   DO I=1,NATOMS-1 
 68:      DO J=I+1, NATOMS 
 69:         IF(LABELS(I,1) /= LABELS(J,1)) THEN 
 70:            IF(K+1<=2*QALCS_NBRHD) THEN 
 71:               NBRHOOD(K) = I; NBRHOOD(K+1) = J 
 72:               K=K+2 
 73:            ELSE 
 74:               WRITE(MYUNIT, '(A)') & 
 75:                    'qalcs_parallel> Neighbourhood size overflow!' 
 76:                
 77:            ENDIF 
 78:         ENDIF 
 79:      ENDDO 
 80:   ENDDO 
 81:   ! 
 82:   IF(K-1 == 2*QALCS_NBRHD) THEN 
 83:      IF(QALCSV) WRITE(MYUNIT,'(A, I6)') & 
 84:           'scan_all_swaps_parallel> Built neighbourhood of size ',& 
 85:           QALCS_NBRHD 
 86:   ELSE 
 87:      WRITE(MYUNIT,'(A, I6, A, I6)') & 
 88:           'scan_all_swaps_parallel> ERROR: neighbourhood size ',& 
 89:           K,' not equal to the expected ', 2*QALCS_NBRHD 
 90:   ENDIF 
 91:   ! 
 92:   COORDSO(1:3*NATOMS,1) = COORDS(1:3*NATOMS,1) 
 93:   LABELSO(1:NATOMS,1) = LABELS(1:NATOMS,1) 
 94:   POTELO = POTEL 
 95:   XMIN(1:3*NATOMS) = COORDS(1:3*NATOMS,1) 
 96:   LMIN(1:NATOMS) = LABELS(1:NATOMS,1) 
 97:   EMIN = POTEL 
 98:   ! 
 99:   IF(QALCSV) WRITE(MYUNIT,'(A)') & 
100:        'scan_all_swaps_parallel> Parallel scan of neighbourhood...' 
101:   ! 
102:   IF(MOD(QALCS_NBRHD,NPAR_GBH)==0) THEN 
103:      NPPN = QALCS_NBRHD/NPAR_GBH 
104:   ELSE 
105:      NPPN = QALCS_NBRHD/NPAR_GBH + 1 ! integer division!!! 
106:   ENDIF 
107:   MYI1 = 2*MYNODE*NPPN+1 
108:   MYI2 = MIN(MYI1-1+2*NPPN, 2*QALCS_NBRHD) 
109:   DO I = MYI1, MYI2, 2 
110:      ! 
111:      !WRITE(MYUNIT,'(A,I6,I6)') & 
112:      !     'scan_all_swaps_parallel> Swapping atoms', & 
113:      !     NBRHOOD(I),NBRHOOD(I+1) 
114:      ! Perform swap and quench 
115:      J=LABELS(NBRHOOD(I),1) 
116:      LABELS(NBRHOOD(I),1) = LABELS(NBRHOOD(I+1),1) 
117:      LABELS(NBRHOOD(I+1),1) = J 
118:      ! 
119:      NQ(1) = NQ(1) + 1 
120:      CALL QUENCH(.FALSE.,1,ITER,TIME,BRUN,QDONE,SCREENC) 
121:      IF(QALCSV) WRITE(MYUNIT,'(A,I10,A,F20.10,A,I5,A,G12.5)') & 
122:           'Qu ',NQ(1),' E= ',POTEL,' steps= ',ITER,' RMS= ',RMS 
123:      ! 
124:      ! Check for lower energy 
125:      IF(POTEL < EMIN - ECONV) THEN 
126:         EMIN = POTEL 
127:         XMIN(:) = COORDS(:,1) 
128:         LMIN(:) = LABELS(:,1) 
129:      ENDIF 
130:      ! 
131:      ! Undo swap and revert state 
132:      !CALL SWAP_LABELS(NBRHOOD(I),NBRHOOD(I+1), 1) 
133:      POTEL = POTELO 
134:      COORDS(:,1) = COORDSO(:,1) 
135:      LABELS(:,1) = LABELSO(:,1) 
136:      ! 
137:   ENDDO 
138:   ! 
139: #ifdef MPI 
140:   ! Gather all parallel quench energies on master node.              
141:   CALL MPI_GATHER(EMIN, 1, MPI_DOUBLE_PRECISION, & 
142:        POTEL_LIST, 1, MPI_DOUBLE_PRECISION, 0, & 
143:        MPI_COMM_WORLD, MPIERR) 
144:   ! 
145:   IF(MYNODE==0) THEN 
146:      ELOWEST=MINVAL(POTEL_LIST) 
147:      IPROC=MINLOC(POTEL_LIST, 1)-1 
148:   ENDIF 
149:   ! 
150:   ! Broadcast IPROC from master 
151:   CALL MPI_BCAST(IPROC,1,MPI_INTEGER,& 
152:        0,MPI_COMM_WORLD,MPIERR) 
153:   ! Broadcast state from IPROC 
154:   CALL MPI_BCAST(XMIN,3*NATOMS,& 
155:        MPI_DOUBLE_PRECISION,IPROC,MPI_COMM_WORLD,MPIERR) 
156:   CALL MPI_BCAST(LMIN,NATOMS,MPI_INTEGER,& 
157:        IPROC,MPI_COMM_WORLD,MPIERR) 
158:   CALL MPI_BCAST(EMIN,1,MPI_DOUBLE_PRECISION,& 
159:        IPROC,MPI_COMM_WORLD,MPIERR) 
160:   IF(TARGET) THEN  
161:      CALL MPI_BCAST(HIT,1,MPI_LOGICAL,& 
162:           IPROC,MPI_COMM_WORLD,MPIERR) 
163:      IF(HIT) WRITE(MYUNIT,'(A,I5)') & 
164:           'QALCS_parallel> Target hit deterministically by node', & 
165:           IPROC+1 
166:   ENDIF 
167: #endif 
168:   ! 
169:   IF(EMIN < POTELO - ECONV) THEN 
170:      COORDS(1:3*NATOMS,1) = XMIN(1:3*NATOMS) 
171:      LABELS(1:NATOMS,1) = LMIN(1:NATOMS) 
172:      POTEL=EMIN 
173:      SWAP=.TRUE. 
174:      IF(QALCSV) WRITE(MYUNIT,'(A, I5)') & 
175:           'scan_all_swaps_parallel> Best improvement on node ',IPROC+1 
176:   ELSE 
177:      SWAP=.FALSE. 
178:      IF(QALCSV) WRITE(MYUNIT,'(A)') & 
179:           'scan_all_swaps_parallel> No improvement found!' 
180:   ENDIF 
181:   ! 
182:   RETURN 
183:   ! 
184: END SUBROUTINE SCAN_ALL_SWAPS_PARALLEL 


r31263/saveit.f 2016-10-12 11:30:27.688833908 +0100 r31262/saveit.f 2016-10-12 11:30:33.668914168 +0100
 93:  93: 
 94:  94: 
 95: ! beh35 > adding stuff to accommodate for mlowest 95: ! beh35 > adding stuff to accommodate for mlowest
 96:  96: 
 97: ! here is how it works normally when you specifiy targets without mlowest modification 97: ! here is how it works normally when you specifiy targets without mlowest modification
 98:        98:       
 99:       IF(TARGET) THEN ! beh35 99:       IF(TARGET) THEN ! beh35
100:       DO J1=1,NTARGETS100:       DO J1=1,NTARGETS
101:          IF (EREAL-TARGETS(J1).LT.ECONV) THEN101:          IF (EREAL-TARGETS(J1).LT.ECONV) THEN
102:             IF (NPAR.LT.2) THEN102:             IF (NPAR.LT.2) THEN
103:                WRITE(MYUNIT,'(2(A,I15),A)') 'saveit> Target hit after ',NQ(NP),' quenches ',NPCALL,' function calls.'103:                WRITE(MYUNIT,'(2(A,I15),A)') 'saveit> Target hit after ',NQTOT,' quenches ',NPCALL,' function calls.'
104:                WRITE(MYUNIT,'(2(A,F20.10))') 'saveit> Energy=',EREAL,' target=',TARGETS(J1)104:                WRITE(MYUNIT,'(2(A,F20.10))') 'saveit> Energy=',EREAL,' target=',TARGETS(J1)
105: !            ELSE105: !            ELSE
106: !               WRITE(MYUNIT,'(A,I0.2,A,I2,2(A,I15),A)') '[',NP,']Target hit in parallel run ',NP,106: !               WRITE(MYUNIT,'(A,I0.2,A,I2,2(A,I15),A)') '[',NP,']Target hit in parallel run ',NP,
107: !     1                    ' after ',NQTOT,' quenches ',NPCALL,' function calls.'107: !     1                    ' after ',NQTOT,' quenches ',NPCALL,' function calls.'
108: !               WRITE(MYUNIT,'(A,I0.2,2(A,F20.10))') '[',NP,']Energy=',EREAL,' target=',TARGETS(J1)108: !               WRITE(MYUNIT,'(A,I0.2,2(A,F20.10))') '[',NP,']Energy=',EREAL,' target=',TARGETS(J1)
109:             ENDIF109:             ENDIF
110:             HIT=.TRUE.110:             HIT=.TRUE.
111:          ENDIF111:          ENDIF
112:       ENDDO112:       ENDDO
113: 113: 


r31263/symmetry.f90 2016-10-12 11:30:28.612846407 +0100 r31262/symmetry.f90 2016-10-12 11:30:33.932917712 +0100
1365:       NQTOT=NQTOT+11365:       NQTOT=NQTOT+1
1366:       NQ(JP)=NQ(JP)+11366:       NQ(JP)=NQ(JP)+1
1367:       CALL QUENCH(.FALSE.,JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)1367:       CALL QUENCH(.FALSE.,JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
1368:       IF (NPAR.GT.1) THEN1368:       IF (NPAR.GT.1) THEN
1369:          WRITE(MYUNIT,'(A,I1,A,I10,A,F20.10,A,I5,A,G12.5,30X,A,F11.1)') &1369:          WRITE(MYUNIT,'(A,I1,A,I10,A,F20.10,A,I5,A,G12.5,30X,A,F11.1)') &
1370:    &                '[',JP,']Qu ',NQ(JP),' E=',POTEL,' steps=',ITERATIONS,' RMS=',RMS,' t=',TIME1370:    &                '[',JP,']Qu ',NQ(JP),' E=',POTEL,' steps=',ITERATIONS,' RMS=',RMS,' t=',TIME
1371:       ELSE1371:       ELSE
1372:          WRITE(MYUNIT,'(A,I10,A,F20.10,A,I5,A,G12.5,30X,A,F11.1)') &1372:          WRITE(MYUNIT,'(A,I10,A,F20.10,A,I5,A,G12.5,30X,A,F11.1)') &
1373:    &                'Qu ',NQ(JP),' E=',POTEL,' steps=',ITERATIONS,' RMS=',RMS,' t=',TIME1373:    &                'Qu ',NQ(JP),' E=',POTEL,' steps=',ITERATIONS,' RMS=',RMS,' t=',TIME
1374:       ENDIF1374:       ENDIF
1375:       !IF (HOMOREFT) THEN <keyword phased out!1375:       IF (HOMOREFT) THEN
1376:       !   CALL HOMOREF(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)1376:          CALL HOMOREF(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
1377:       !ENDIF1377:       ENDIF
1378:    ! Save the best structure1378:    ! Save the best structure
1379:       IF (LEBEST-POTEL.GT.ECONV) THEN1379:       IF (LEBEST-POTEL.GT.ECONV) THEN
1380:          LEBEST=POTEL1380:          LEBEST=POTEL
1381:          QBEST(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)1381:          QBEST(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)
1382:          VATBEST(1:NATOMS)=VAT(1:NATOMS,JP)1382:          VATBEST(1:NATOMS)=VAT(1:NATOMS,JP)
1383:          QBCHANGED=.TRUE.1383:          QBCHANGED=.TRUE.
1384:          QBORDERED=.TRUE.1384:          QBORDERED=.TRUE.
1385:          QBCORE=NCOREREAL1385:          QBCORE=NCOREREAL
1386: !        WRITE(MYUNIT,'(A,2L5,I6)') 'symmetry> D QBORDERED,QBCHANGED,QBCORE=',QBORDERED,QBCHANGED,QBCORE1386: !        WRITE(MYUNIT,'(A,2L5,I6)') 'symmetry> D QBORDERED,QBCHANGED,QBCORE=',QBORDERED,QBCHANGED,QBCORE
1387: 1387: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0