hdiff output

r32062/finalio.f90 2017-03-10 15:30:13.677252368 +0000 r32061/finalio.f90 2017-03-10 15:30:14.461262735 +0000
 24:   USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_FINISH, AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, & 24:   USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_FINISH, AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, &
 25:        AMBER12_WRITE_XYZ 25:        AMBER12_WRITE_XYZ
 26:   USE OPEP_INTERFACE_MOD, ONLY: OPEP_FINISH, OPEP_WRITE_PDB 26:   USE OPEP_INTERFACE_MOD, ONLY: OPEP_FINISH, OPEP_WRITE_PDB
 27:   USE QMODULE 27:   USE QMODULE
 28:   USE MODCHARMM 28:   USE MODCHARMM
 29:   USE AMHGLOBALS, ONLY:NMRES,IRES 29:   USE AMHGLOBALS, ONLY:NMRES,IRES
 30:   USE BGUPMOD 30:   USE BGUPMOD
 31:   USE PERMU 31:   USE PERMU
 32:   USE MODHESS, ONLY : HESS, MASSWT 32:   USE MODHESS, ONLY : HESS, MASSWT
 33:   USE CONVEX_POLYHEDRA_MODULE, ONLY: VIEW_POLYHEDRA 33:   USE CONVEX_POLYHEDRA_MODULE, ONLY: VIEW_POLYHEDRA
 34:   USE LJ_GAUSS_MOD, ONLY: VIEW_LJ_GAUSS 34:   
 35:  
 36:   IMPLICIT NONE 35:   IMPLICIT NONE
 37:  36:   
 38:   !   MCP 37:   !   MCP
 39:   INTEGER III, I3,  GLY_COUNT, ID, NUMCRD, NUMPRO, NCPHST, GETUNIT, AMHUNIT1, LUNIT, NSYMOPS, NCLOSE 38:   INTEGER III, I3,  GLY_COUNT, ID, NUMCRD, NUMPRO, NCPHST, GETUNIT, AMHUNIT1, AMHUNIT2, LUNIT, NSYMOPS, NCLOSE
 40:   INTEGER J1, J2, J3, J4, J5, MYUNIT2, I1, NDUMMY, MYUNIT3, NRBS1, NRBS2, LJGOUNIT, MJ2, REORDERUNIT 39:   INTEGER J1, J2, J3, J4, J5, MYUNIT2, I1, NDUMMY, MYUNIT3, NC, NRBS1, NRBS2, MYUNIT4, LJGOUNIT, MJ2, REORDERUNIT
 41:   DOUBLE PRECISION RBCOORDS(NRBSITES*3), DCOORDS(3*NATOMS), EDUMMY, ITDET, DIST 40:   DOUBLE PRECISION RBCOORDS(NRBSITES*3), DCOORDS(3*NATOMS), EDUMMY, ITDET, DIST
 42:   ! 41:   !
 43:   !ds656> Symmetry detection for clusters on a substrate... 42:   !ds656> Symmetry detection for clusters on a substrate...
 44:   INTEGER :: ISITES(MIEF_NSITES), NSITES, COMBO_N 43:   INTEGER :: ISITES(MIEF_NSITES), NSITES, COMBO_N
 45:   DOUBLE PRECISION :: DISTS(MIEF_NSITES) 44:   DOUBLE PRECISION :: SITES(3*MIEF_NSITES), DISTS(MIEF_NSITES)
 46:   INTEGER, ALLOCATABLE :: COMBO_LABELS(:) 45:   INTEGER, ALLOCATABLE :: COMBO_LABELS(:)
 47:   DOUBLE PRECISION, ALLOCATABLE :: COMBO_COORDS(:) 46:   DOUBLE PRECISION, ALLOCATABLE :: COMBO_COORDS(:)
 48:   !<ds656 47:   !<ds656
 49:     ! 48:     !
 50:   DOUBLE PRECISION P(3), DU(3), RMI(3,3), DRMI(3,3), PI, PHI, THT 49:   DOUBLE PRECISION P3(3,3), P(3), DU(3), RMI(3,3), DRMI(3,3), PI, PHI, THT, CHI
 51:   !DOUBLE PRECISION, ALLOCATABLE :: XCOORDS(:), YCOORDS(:) 50:   !DOUBLE PRECISION, ALLOCATABLE :: XCOORDS(:), YCOORDS(:)
 52:   CHARACTER(LEN=13) J1CHAR,J1CHAR2                  !for gay-berne output files 51:   CHARACTER(LEN=13) J1CHAR,J1CHAR2                  !for gay-berne output files
 53:   CHARACTER(LEN=6) ZSTR 52:   CHARACTER(LEN=6) ZSTR
 54: !  CHARACTER(LEN=2)  DUMMYSTR 53:   CHARACTER(LEN=2)  DUMMYSTR
 55:   DOUBLE PRECISION EulerPhi,EulerPsi,EulerTheta,EulerThetadeg,EulerPhiDeg,EulerPsiDeg  ! Euler angles for ellipsoids of revolution 54:   DOUBLE PRECISION EulerPhi,EulerPsi,EulerTheta,EulerThetadeg,EulerPhiDeg,EulerPsiDeg  ! Euler angles for ellipsoids of revolution
 56:  55:   
 57:   DOUBLE PRECISION EPS2, RAD, HEIGHT,sumx,sumy,sumz,CM(3),SYMOPS(120,3,3) 56:   DOUBLE PRECISION EPS2, RAD, HEIGHT,sumx,sumy,sumz,CM(3),SYMOPS(120,3,3)
 58:   LOGICAL :: GTEST, ADDDONE(NATOMS) 57:   LOGICAL :: GTEST, ADDDONE(NATOMS), ADDORDER(NATOMS)
 59:   COMMON /CAPS/ EPS2, RAD, HEIGHT 58:   COMMON /CAPS/ EPS2, RAD, HEIGHT
 60:  59:   
 61:   CHARACTER(LEN=4) :: FPGRP 60:   CHARACTER(LEN=4) :: FPGRP
 62:   CHARACTER(LEN=6) :: CRMS 61:   CHARACTER(LEN=6) :: CRMS
 63:   CHARACTER(LEN=20) :: MYFILENAME2, ISTR, DBNUM!, MYFILENAME3 62:   CHARACTER(LEN=20) :: MYFILENAME2, ISTR, DBNUM, MYFILENAME3
 64:   CHARACTER(LEN=15), ALLOCATABLE :: DBNAME(:) 63:   CHARACTER(LEN=15), ALLOCATABLE :: DBNAME(:)
 65:   CHARACTER(LEN=80) :: tempstring 64:   CHARACTER(LEN=80) :: tempstring
 66:  65:   
 67:   CHARACTER(LEN=6) :: J1_STRING, MYNODE_STRING 66:   CHARACTER(LEN=6) :: J1_STRING, MYNODE_STRING
 68:  67:   
 69:   !  AMH 68:   !  AMH
 70:   CHARACTER(LEN=3) :: RES_TYPE 69:   CHARACTER(LEN=3) :: RES_TYPE
 71:   CHARACTER(LEN=2) :: ATOM_TYPE 70:   CHARACTER(LEN=2) :: ATOM_TYPE
 72:   CHARACTER*1 COUNTTT 71:   CHARACTER*1 COUNTTT
 73:   INTEGER COUNTT 72:   INTEGER COUNTT
 74:   DOUBLE PRECISION  PPPCORD(NMRES*3*3,3,3,5) 73:   DOUBLE PRECISION  PPPCORD(NMRES*3*3,3,3,5)
 75:   EXTERNAL NUM_TO_CHAR 74:   EXTERNAL NUM_TO_CHAR
 76:   DOUBLE PRECISION TEND, DMIN 75:   DOUBLE PRECISION TEND, DMIN
 77: !  INTEGER QDONE, NP, ITERATIONS, CLOSEST(NADDTARGET), BRUN, MYUNIT4 76:   INTEGER ITERATIONS, BRUN,QDONE, NP, CLOSEST(NADDTARGET)
 78: !  DOUBLE PRECISION SCREENC(3*NATOMSALLOC), TIME 77:   DOUBLE PRECISION SCREENC(3*NATOMSALLOC), TIME
 79:   CHARACTER(LEN=2) :: ALPHABET(26)=(/'ja','jb','jc','jd','je','jf','jg','jh','ji','jj','jk','jl', & 78:   CHARACTER(LEN=2) :: ALPHABET(26)=(/'ja','jb','jc','jd','je','jf','jg','jh','ji','jj','jk','jl', &
 80:   &                                  'jm','jn','jo','jp','jq','jr','js','jt','ju','jv','jw','jx','jy','jz'/) 79:   &                                  'jm','jn','jo','jp','jq','jr','js','jt','ju','jv','jw','jx','jy','jz'/)
 81:  80:   
 82:   ! ds656> Extra stuff for HSA in OPTIM's min.data format 81:   ! ds656> Extra stuff for HSA in OPTIM's min.data format
 83:   INTEGER :: NUM_ZERO_EVS, INFO 82:   INTEGER :: NUM_ZERO_EVS, INFO
 84:   DOUBLE PRECISION :: LOG_PROD, EVALS(3*NATOMS), WORK(9*NATOMS) 83:   DOUBLE PRECISION :: LOG_PROD, EVALS(3*NATOMS), WORK(9*NATOMS)
 85:  84:   
 86:  85:   
 87:   PI = 4.D0*DATAN(1.D0) 86:   PI = 4.D0*DATAN(1.D0)
 88:  87:   
 89:   !ds656> test 88:   !ds656> test
 90:   !write(*,*) "finalio> START, PI=", PI 89:   !write(*,*) "finalio> START, PI=", PI
 91:  90:   
 92:   IF (PERMOPT) THEN 91:   IF (PERMOPT) THEN
 93:      NSAVE=2 92:      NSAVE=2
 94:      QMIN(2)=0.0D0 93:      QMIN(2)=0.0D0
 95:      QMINP(2,:)=FIN(:) 94:      QMINP(2,:)=FIN(:)
 96:      FF(:)=0 95:      FF(:)=0
 97:   ENDIF 96:   ENDIF
 98:  97:   
 99:   NUMPRO = 1 98:   NUMPRO = 1
100:   NUMCRD = 3 99:   NUMCRD = 3
101: 100:   
102:   ALLOCATE(DBNAME(NSAVE))101:   ALLOCATE(DBNAME(NSAVE))
103: 102:   
104:   IF (DEBUG) WRITE(MYUNIT,'(A,3I6)') ' in finalio MYNODE,MYUNIT,NSAVE=',MYNODE,MYUNIT,NSAVE !jdf43>103:   IF (DEBUG) WRITE(MYUNIT,'(A,3I6)') ' in finalio MYNODE,MYUNIT,NSAVE=',MYNODE,MYUNIT,NSAVE !jdf43>
105:   DO J1=1,NSAVE104:   DO J1=1,NSAVE
106:      WRITE(DBNUM,*) J1105:      WRITE(DBNUM,*) J1
107:      DBNAME(J1)='dbase.'//TRIM(ADJUSTL(DBNUM))106:      DBNAME(J1)='dbase.'//TRIM(ADJUSTL(DBNUM))
108:   ENDDO107:   ENDDO
109:   DCOORDS(1:3*NATOMS)=0.0D0108:   DCOORDS(1:3*NATOMS)=0.0D0
110:   !      IF (AMH) THEN109:   !      IF (AMH) THEN
111:   !         CALL WALESAMH_FINALIO110:   !         CALL WALESAMH_FINALIO
112:   !         STOP111:   !         STOP
113:   !      ENDIF112:   !      ENDIF
114: 113:   
115:   IF (AMHT) THEN114:   IF (AMHT) THEN
116:      AMHUNIT1=GETUNIT()115:      AMHUNIT1=GETUNIT()
117:      OPEN(UNIT=AMHUNIT1,FILE='movie_gmin',STATUS='UNKNOWN')116:      OPEN(UNIT=AMHUNIT1,FILE='movie_gmin',STATUS='UNKNOWN')
118:      WRITE(AMHUNIT1,334)NMRES,NUMCRD,NUMPRO,NSAVE117:      WRITE(AMHUNIT1,334)NMRES,NUMCRD,NUMPRO,NSAVE
119: 334  FORMAT(4(I8,1X),' NMRES NMCRD NUMPRO NMSNAP')118: 334  FORMAT(4(I8,1X),' NMRES NMCRD NUMPRO NMSNAP')
120:   ENDIF119:   ENDIF
121: 120:   
122:   IF (DMACRYST) THEN121:   IF (DMACRYST) THEN
123:      CALL DMACRYS_DUMP_LOWEST122:      CALL DMACRYS_DUMP_LOWEST
124:   ENDIF123:   ENDIF
125: 124:   
126:   IF(USERPOTT) THEN125:   IF(USERPOTT) THEN
127:      CALL USERPOT_DUMP_LOWEST126:      CALL USERPOT_DUMP_LOWEST
128:   ENDIF127:   ENDIF
129: 128:   
130: 129:   
131:   IF (MPIT.OR.GBHT) THEN130:   IF (MPIT.OR.GBHT) THEN
132:      WRITE (ISTR, '(I10)') MYNODE+1131:      WRITE (ISTR, '(I10)') MYNODE+1
133:      MYUNIT2=GETUNIT()132:      MYUNIT2=GETUNIT()
134:      MYFILENAME2="lowest."//TRIM(ADJUSTL(ISTR))133:      MYFILENAME2="lowest."//TRIM(ADJUSTL(ISTR))
135:      !        WRITE(MYUNIT,'(A,I6,A)') ' in finalio MYUNIT2,MYFILENAME2=',MYUNIT2,TRIM(ADJUSTL(MYFILENAME2))134:      !        WRITE(MYUNIT,'(A,I6,A)') ' in finalio MYUNIT2,MYFILENAME2=',MYUNIT2,TRIM(ADJUSTL(MYFILENAME2))
136:      OPEN(MYUNIT2,FILE=TRIM(ADJUSTL(MYFILENAME2)), STATUS="UNKNOWN", FORM="FORMATTED")135:      OPEN(MYUNIT2,FILE=TRIM(ADJUSTL(MYFILENAME2)), STATUS="UNKNOWN", FORM="FORMATTED")
137:      IF (CHRMMT) THEN136:      IF (CHRMMT) THEN
138:         DO J1=1,NSAVE137:         DO J1=1,NSAVE
139:            WRITE(DBNUM,*) J1138:            WRITE(DBNUM,*) J1
140:            DBNAME(J1)='dbase.'//TRIM(ADJUSTL(ISTR))//'.'//TRIM(ADJUSTL(DBNUM))139:            DBNAME(J1)='dbase.'//TRIM(ADJUSTL(ISTR))//'.'//TRIM(ADJUSTL(DBNUM))
141:         ENDDO140:         ENDDO
142:      ENDIF141:      ENDIF
143:      IF (CSMT.AND.(.NOT.SYMMETRIZECSM)) THEN142:      IF (CSMT.AND.(.NOT.SYMMETRIZECSM)) THEN
144:         MYUNIT3=GETUNIT()143:         MYUNIT3=GETUNIT()
145: !        MYFILENAME3="CSMav."//trim(adjustl(istr))144:         MYFILENAME3="CSMav."//trim(adjustl(istr))
146:      ENDIF145:      ENDIF
147:   ELSE146:   ELSE
148:      MYUNIT2=GETUNIT()147:      MYUNIT2=GETUNIT()
149:      ! hk286148:      ! hk286
150:      IF (RIGIDINIT) THEN149:      IF (RIGIDINIT) THEN
151:         OPEN(MYUNIT2,FILE='GRlowest',STATUS='UNKNOWN')150:         OPEN(MYUNIT2,FILE='GRlowest',STATUS='UNKNOWN')
152:      ELSE151:      ELSE
153:         OPEN(MYUNIT2,FILE='lowest',STATUS='UNKNOWN')152:         OPEN(MYUNIT2,FILE='lowest',STATUS='UNKNOWN')
154:      ENDIF153:      ENDIF
155:      ! hk286154:      ! hk286
156:      IF (CSMT.AND.(.NOT.SYMMETRIZECSM)) THEN155:      IF (CSMT.AND.(.NOT.SYMMETRIZECSM)) THEN
157:         MYUNIT3=26156:         MYUNIT3=26 
158:         OPEN(MYUNIT3,FILE='CSMav.xyz',STATUS='UNKNOWN')157:         OPEN(MYUNIT3,FILE='CSMav.xyz',STATUS='UNKNOWN')
159:      ENDIF158:      ENDIF
160:   ENDIF159:   ENDIF
161: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!160: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
162: !161: !
163: ! debug DJW162: ! debug DJW
164: !163: !
165: !        MYUNIT3=GETUNIT()164: !        MYUNIT3=GETUNIT()
166: !        OPEN(MYUNIT3,FILE='GCBHlowest',STATUS='UNKNOWN')165: !        OPEN(MYUNIT3,FILE='GCBHlowest',STATUS='UNKNOWN')
167: !        MYUNIT4=GETUNIT()166: !        MYUNIT4=GETUNIT()
168: !        OPEN(MYUNIT4,FILE='GCBHlowest.new',STATUS='UNKNOWN')167: !        OPEN(MYUNIT4,FILE='GCBHlowest.new',STATUS='UNKNOWN')
169: ! 777    READ(MYUNIT3,*,END=888) J1168: ! 777    READ(MYUNIT3,*,END=888) J1
170: !        NATOMS=J1169: !        NATOMS=J1
171: !        READ(MYUNIT3,*) DUMMYSTR170: !        READ(MYUNIT3,*) DUMMYSTR
172: !        NP=1171: !        NP=1
173: !        DO J2=1,J1172: !        DO J2=1,J1
174: !           READ(MYUNIT3,'(A2,1X,3G20.10)') DUMMYSTR,COORDS(3*(J2-1)+1:3*(J2-1)+3,NP)173: !           READ(MYUNIT3,'(A2,1X,3G20.10)') DUMMYSTR,COORDS(3*(J2-1)+1:3*(J2-1)+3,NP)
175: !        ENDDO174: !        ENDDO
176: !175: ! 
177: !        CALL QUENCH(.TRUE.,NP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)176: !        CALL QUENCH(.TRUE.,NP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
178: !        WRITE(MYUNIT4,'(I10)') J1177: !        WRITE(MYUNIT4,'(I10)') J1
179: !        WRITE(MYUNIT4,'(A,2G20.10)') 'energy=',QPE(J1),QENERGIES(J1)178: !        WRITE(MYUNIT4,'(A,2G20.10)') 'energy=',QPE(J1),QENERGIES(J1)
180: !        DO J2=1,J1179: !        DO J2=1,J1
181: !           WRITE(MYUNIT4,'(A2,1X,3G20.10)') 'AX',QCOORDINATES(J1,3*(J2-1)+1:3*(J2-1)+3)180: !           WRITE(MYUNIT4,'(A2,1X,3G20.10)') 'AX',QCOORDINATES(J1,3*(J2-1)+1:3*(J2-1)+3)
182: !        ENDDO181: !        ENDDO
183: !182: ! 
184: !        GOTO 777183: !        GOTO 777
185: ! 888    CLOSE(MYUNIT3)184: ! 888    CLOSE(MYUNIT3)
186: !        CLOSE(MYUNIT4)185: !        CLOSE(MYUNIT4)
187: !        STOP186: !        STOP
188: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!187: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
189:   IF (GCBHT) THEN188:   IF (GCBHT) THEN
190:      MYUNIT3=GETUNIT()189:      MYUNIT3=GETUNIT()
191:      OPEN(MYUNIT3,FILE='GCBHlowest',STATUS='UNKNOWN')190:      OPEN(MYUNIT3,FILE='GCBHlowest',STATUS='UNKNOWN')
192:      DO J1=1,GCNATOMS191:      DO J1=1,GCNATOMS
193:         IF (QENERGIES(J1).LT.HUGE(1.0D0)/10.0) THEN192:         IF (QENERGIES(J1).LT.HUGE(1.0D0)/10.0) THEN
194:            WRITE(MYUNIT3,'(I10)') J1193:            WRITE(MYUNIT3,'(I10)') J1
195:            WRITE(MYUNIT3,'(A,2G20.10)') 'energy=',QPE(J1),QENERGIES(J1)194:            WRITE(MYUNIT3,'(A,2G20.10)') 'energy=',QPE(J1),QENERGIES(J1)
196:            DO J2=1,J1195:            DO J2=1,J1
197:               WRITE(MYUNIT3,'(A2,1X,3G20.10)') 'AX',QCOORDINATES(J1,3*(J2-1)+1:3*(J2-1)+3)196:               WRITE(MYUNIT3,'(A2,1X,3G20.10)') 'AX',QCOORDINATES(J1,3*(J2-1)+1:3*(J2-1)+3)
198:            ENDDO197:            ENDDO
199:         ENDIF198:         ENDIF
200:      ENDDO199:      ENDDO
201:      CLOSE(MYUNIT3)200:      CLOSE(MYUNIT3)
202:   ENDIF201:   ENDIF
203: 202: 
204:   IF (LJADD3T) THEN203:   IF (LJADD3T) THEN 
205:      LJGOUNIT=GETUNIT()204:      LJGOUNIT=GETUNIT()
206:      OPEN(LJGOUNIT,FILE='ljadd3.xyz',STATUS='UNKNOWN')205:      OPEN(LJGOUNIT,FILE='ljadd3.xyz',STATUS='UNKNOWN')
207:   ENDIF206:   ENDIF
208:   IF (REORDERADDT) THEN207:   IF (REORDERADDT) THEN 
209:      REORDERUNIT=GETUNIT()208:      REORDERUNIT=GETUNIT()
210:      OPEN(REORDERUNIT,FILE='reorder.xyz',STATUS='UNKNOWN')209:      OPEN(REORDERUNIT,FILE='reorder.xyz',STATUS='UNKNOWN')
211:   ENDIF210:   ENDIF
212: 211:   
213:   savemin: DO J1=1,NSAVE212:   savemin: DO J1=1,NSAVE
214:      !213:      !
215:      NATOMS=QMINNATOMS(J1)214:      NATOMS=QMINNATOMS(J1)
216:      IF (AMHT) THEN215:      IF (AMHT) THEN
217:         COUNTT=J1216:         COUNTT=J1
218:         CALL NUM_TO_CHAR(COUNTT,COUNTTT)217:         CALL NUM_TO_CHAR(COUNTT,COUNTTT)
219:         OPEN(UNIT=27,FILE='movie_gmin.'//COUNTTT//'.pdb',STATUS='UNKNOWN')218:         OPEN(UNIT=27,FILE='movie_gmin.'//COUNTTT//'.pdb',STATUS='UNKNOWN')
220:      ENDIF219:      ENDIF
221: 220:      
222:      IF (RGCL2.OR.ARNO) THEN221:      IF (RGCL2.OR.ARNO) THEN
223:         WRITE(MYUNIT2,*) NATOMS+2222:         WRITE(MYUNIT2,*) NATOMS+2
224:      ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT) THEN223:      ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT) THEN
225:         WRITE(MYUNIT2,*) NATOMS/2224:         WRITE(MYUNIT2,*) NATOMS/2
226:      ELSE IF (AMHT) THEN225:      ELSE IF (AMHT) THEN
227:         WRITE(MYUNIT2,*) NMRES*3226:         WRITE(MYUNIT2,*) NMRES*3
228:      ELSE IF (.NOT. PRINT_MINDATA) THEN227:      ELSE IF (.NOT. PRINT_MINDATA) THEN
229:         WRITE(MYUNIT2,*) NATOMS228:         WRITE(MYUNIT2,*) NATOMS
230:         !ds656> With PRINT_MINDATA we want to discard "minima"229:         !ds656> With PRINT_MINDATA we want to discard "minima"
231:         ! with negative eigenvalues, and so we print the number230:         ! with negative eigenvalues, and so we print the number
257:               ISITES(J4) = J4256:               ISITES(J4) = J4
258:               DISTS(J4) = 0.0257:               DISTS(J4) = 0.0
259:               DO J5=1,3258:               DO J5=1,3
260:                  DISTS(J4) = DISTS(J4) + &259:                  DISTS(J4) = DISTS(J4) + &
261:                       (MIEF_SITES(J4,J5)-CM(J5))**2260:                       (MIEF_SITES(J4,J5)-CM(J5))**2
262:               ENDDO261:               ENDDO
263:               DISTS(J4) = DSQRT(DISTS(J4))262:               DISTS(J4) = DSQRT(DISTS(J4))
264:               J5=J4263:               J5=J4
265:               sort_sites: DO WHILE(J5>1)264:               sort_sites: DO WHILE(J5>1)
266:                  IF(DISTS(J5) < DISTS(J5-1)) THEN265:                  IF(DISTS(J5) < DISTS(J5-1)) THEN
267:                     NDUMMY = ISITES(J5)266:                     NDUMMY = ISITES(J5) 
268:                     ISITES(J5) = ISITES(J5-1)267:                     ISITES(J5) = ISITES(J5-1)
269:                     ISITES(J5-1) = NDUMMY268:                     ISITES(J5-1) = NDUMMY
270:                     EDUMMY = DISTS(J5)269:                     EDUMMY = DISTS(J5)
271:                     DISTS(J5) = DISTS(J5-1)270:                     DISTS(J5) = DISTS(J5-1)
272:                     DISTS(J5-1) = EDUMMY271:                     DISTS(J5-1) = EDUMMY
273:                     J5=J5-1272:                     J5=J5-1
274:                  ELSE273:                  ELSE
275:                     EXIT sort_sites274:                     EXIT sort_sites
276:                  ENDIF275:                  ENDIF
277:               ENDDO sort_sites276:               ENDDO sort_sites
278:            ENDDO277:            ENDDO
279:            !278:            ! 
280:            !ds656> Extract a small subset of field sites.279:            !ds656> Extract a small subset of field sites. 
281:            ! While coding this I had graphene-like lattice in280:            ! While coding this I had graphene-like lattice in
282:            ! in mind, in which case the fragment can be small.281:            ! in mind, in which case the fragment can be small.
283:            ! Complications may arise for other substrates...282:            ! Complications may arise for other substrates...
284:            J5=3283:            J5=3 
285:            count_sites: DO J4=4,MIEF_NSITES284:            count_sites: DO J4=4,MIEF_NSITES
286:               IF(ABS(DISTS(J4)-DISTS(J4-1)) < PGSYMTOLS(2)) THEN285:               IF(ABS(DISTS(J4)-DISTS(J4-1)) < PGSYMTOLS(2)) THEN
287:                  J5=J5+1286:                  J5=J5+1
288:               ELSE287:               ELSE
289:                  EXIT count_sites288:                  EXIT count_sites
290:               ENDIF289:               ENDIF
291:            ENDDO count_sites290:            ENDDO count_sites
292:            NSITES=J5291:            NSITES=J5
293:            !292:            !
294:            COMBO_N=NATOMS+NSITES293:            COMBO_N=NATOMS+NSITES
309:            COMBO_LABELS(NSITES+1:COMBO_N) = QMINT(J1,1:NATOMS)308:            COMBO_LABELS(NSITES+1:COMBO_N) = QMINT(J1,1:NATOMS)
310:            COMBO_COORDS(3*NSITES+1:3*COMBO_N) = &309:            COMBO_COORDS(3*NSITES+1:3*COMBO_N) = &
311:                 QMINP(J1,1:3*NATOMS)310:                 QMINP(J1,1:3*NATOMS)
312:            !311:            !
313:            WRITE(MYUNIT,'(A,I4,A)') &312:            WRITE(MYUNIT,'(A,I4,A)') &
314:                 'finalio> Diagnosing point group including ', &313:                 'finalio> Diagnosing point group including ', &
315:                 NSITES,' MIEF sites.'314:                 NSITES,' MIEF sites.'
316:            !315:            !
317:            CALL PGSYM( COMBO_N, COMBO_COORDS(1:3*COMBO_N), &316:            CALL PGSYM( COMBO_N, COMBO_COORDS(1:3*COMBO_N), &
318:                 COMBO_LABELS(1:COMBO_N), PGSYMTOLS, 2, J3,&317:                 COMBO_LABELS(1:COMBO_N), PGSYMTOLS, 2, J3,&
319:                 CM, NSYMOPS, SYMOPS, FPGRP )318:                 CM, NSYMOPS, SYMOPS, FPGRP )              
320:            !319:            !
321:            DEALLOCATE(COMBO_COORDS, COMBO_LABELS)320:            DEALLOCATE(COMBO_COORDS, COMBO_LABELS)
322:            !321:            !
323:         ENDIF322:         ENDIF
324:         !323:         !
325:         IF(PRINT_PTGRP) THEN ! Just append point group324:         IF(PRINT_PTGRP) THEN ! Just append point group
326:            WRITE(MYUNIT2,99) J1, QMIN(J1), FF(J1), &325:            WRITE(MYUNIT2,99) J1, QMIN(J1), FF(J1), &
327:                 NPCALL_QMIN(J1), FPGRP, NSYMOPS326:                 NPCALL_QMIN(J1), FPGRP, NSYMOPS
328: 99         FORMAT('Energy of minimum ',I6,'=',G20.10, &327: 99         FORMAT('Energy of minimum ',I6,'=',G20.10, &
329:                 ' first found at step ',I8,' after ', &328:                 ' first found at step ',I8,' after ', &
340:            IF(NSPECIES(0)>1) THEN ! ds656> this also sets ATMASS339:            IF(NSPECIES(0)>1) THEN ! ds656> this also sets ATMASS
341:               CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)340:               CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)
342:               CALL SET_LABELS(QMINT(J1,1:NATOMS),1)341:               CALL SET_LABELS(QMINT(J1,1:NATOMS),1)
343:            ENDIF342:            ENDIF
344:            CALL POTENTIAL(QMINP(J1,1:3*NATOMS), EVALS, EDUMMY, &343:            CALL POTENTIAL(QMINP(J1,1:3*NATOMS), EVALS, EDUMMY, &
345:                 .TRUE., .TRUE.) ! EVALS is dummy for gradient here344:                 .TRUE., .TRUE.) ! EVALS is dummy for gradient here
346:            IF(ABS(EDUMMY-QMIN(J1)) > ECONV) THEN345:            IF(ABS(EDUMMY-QMIN(J1)) > ECONV) THEN
347:               WRITE(MYUNIT,'(A)') &346:               WRITE(MYUNIT,'(A)') &
348:                    'finalio> WARNING: Unexpected energy value! Configuration changed on final requench?'347:                    'finalio> WARNING: Unexpected energy value! Configuration changed on final requench?'
349:            ENDIF348:            ENDIF
350:            CALL MASSWT() ! Weight the Hessian by mass.349:            CALL MASSWT() ! Weight the Hessian by mass. 
351:            !350:            !
352:            WRITE(MYUNIT, '(A)') &351:            WRITE(MYUNIT, '(A)') &
353:                 'finalio> Using DSYEV to calculate Hessian eigenvalues'352:                 'finalio> Using DSYEV to calculate Hessian eigenvalues'
354:            CALL DSYEV('N','U',3*NATOMS,HESS,3*NATOMS,EVALS,&353:            CALL DSYEV('N','U',3*NATOMS,HESS,3*NATOMS,EVALS,&
355:                 WORK,9*NATOMS,INFO)354:                 WORK,9*NATOMS,INFO)
356:            IF(INFO /= 0) THEN355:            IF(INFO /= 0) THEN
357:               WRITE(MYUNIT,'(A,I2)') &356:               WRITE(MYUNIT,'(A,I2)') &
358:                    'finalio> ***WARNING*** DSYEV error status', INFO357:                    'finalio> ***WARNING*** DSYEV error status', INFO
359:            ENDIF358:            ENDIF
360:            J3=0359:            J3=0
361:            WRITE(MYUNIT,'(A,I6,A)',ADVANCE='NO') &360:            WRITE(MYUNIT,'(A,I6,A)',ADVANCE='NO') &
362:                 'finalio> ',3*NATOMS,' Hessian eigenvalues:'361:                 'finalio> ',3*NATOMS,' Hessian eigenvalues:'
363:            DO J2=1,3*NATOMS362:            DO J2=1,3*NATOMS
364:               WRITE(MYUNIT,'(1X, E12.5E2)',ADVANCE='NO') EVALS(J2)363:               WRITE(MYUNIT,'(1X, E12.5E2)',ADVANCE='NO') EVALS(J2)
365:               IF(EVALS(J2) < - ABS(HESS_EIGEN_TOL)) THEN364:               IF(EVALS(J2) < - ABS(HESS_EIGEN_TOL)) THEN
366:                  WRITE(MYUNIT,*)365:                  WRITE(MYUNIT,*)
367:                  WRITE(MYUNIT,'(A)') &366:                  WRITE(MYUNIT,'(A)') &
368:                       'finalio> WARNING: -ve eigenvalue => not writing to ''lowest''!'367:                       'finalio> WARNING: -ve eigenvalue => not writing to ''lowest''!'
369:                  CYCLE savemin368:                  CYCLE savemin
370:               ELSEIF(ABS(EVALS(J2)) < ABS(HESS_EIGEN_TOL)) THEN369:               ELSEIF(ABS(EVALS(J2)) < ABS(HESS_EIGEN_TOL)) THEN 
371:                  J3=J3+1370:                  J3=J3+1
372:               ENDIF371:               ENDIF
373:            ENDDO372:            ENDDO
374:            WRITE(MYUNIT,*)373:            WRITE(MYUNIT,*)
375:            IF(J3 /= NUM_ZERO_EVS) THEN374:            IF(J3 /= NUM_ZERO_EVS) THEN
376:               WRITE(MYUNIT,'(A,I6,A)') &375:               WRITE(MYUNIT,'(A,I6,A)') &
377:                    'finalio> WARNING: Unexpected number of zero eigenvalues: ',J3,' => not writing to ''lowest''!'376:                    'finalio> WARNING: Unexpected number of zero eigenvalues: ',J3,' => not writing to ''lowest''!'
378:            ELSE377:            ELSE
379:               WRITE(MYUNIT,'(A,I6,A)') &378:               WRITE(MYUNIT,'(A,I6,A)') &
380:                    'finalio> Detected ', J3,' zero eigenvalues.'379:                    'finalio> Detected ', J3,' zero eigenvalues.'
381:            ENDIF380:            ENDIF
382:            ! Log prod. of frequencies is log prod. of sqrt of evals381:            ! Log prod. of frequencies is log prod. of sqrt of evals 
383:            LOG_PROD = &382:            LOG_PROD = &
384:                 SUM(DLOG(EVALS(NUM_ZERO_EVS+1:3*NATOMS)))383:                 SUM(DLOG(EVALS(NUM_ZERO_EVS+1:3*NATOMS)))
385:            WRITE(MYUNIT,'(A,I6,A,G20.10)') &384:            WRITE(MYUNIT,'(A,I6,A,G20.10)') &
386:                 'finalio> Log product of ', &385:                 'finalio> Log product of ', &
387:                 3*NATOMS-NUM_ZERO_EVS,&386:                 3*NATOMS-NUM_ZERO_EVS,&
388:                 ' non-zero eigenvalues: ', &387:                 ' non-zero eigenvalues: ', &
389:                 LOG_PROD388:                 LOG_PROD
390:            CALL MYINERTIA(QMINP(J1,1:3*NATOMS),ITDET)389:            CALL MYINERTIA(QMINP(J1,1:3*NATOMS),ITDET)
391:            WRITE(MYUNIT,'(A,G20.10)') &390:            WRITE(MYUNIT,'(A,G20.10)') &
392:                 'finalio> Determinant of inertia tensor: ', ITDET391:                 'finalio> Determinant of inertia tensor: ', ITDET
412:         WRITE(MYUNIT2,'(A,F20.10)') 'Cl 0.0 0.0 ',-0.995D0411:         WRITE(MYUNIT2,'(A,F20.10)') 'Cl 0.0 0.0 ',-0.995D0
413:         WRITE(MYUNIT2,60) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))412:         WRITE(MYUNIT2,60) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))
414: 60      FORMAT('AR ',3F20.10)413: 60      FORMAT('AR ',3F20.10)
415:      ELSE IF (AMHT) THEN414:      ELSE IF (AMHT) THEN
416:         !415:         !
417:         !   OUTPUT CORDS FOR LOWEST IN X,Y,Z FORMAT416:         !   OUTPUT CORDS FOR LOWEST IN X,Y,Z FORMAT
418:         !   OUTPUT CORDS FOR MOVIE_GMIN IN MOVIESEG FORMAT417:         !   OUTPUT CORDS FOR MOVIE_GMIN IN MOVIESEG FORMAT
419:         !418:         !
420:         WRITE(AMHUNIT1,683)NUMCRD,j1,NUMCRD,REAL(NUMCRD),NUMCRD419:         WRITE(AMHUNIT1,683)NUMCRD,j1,NUMCRD,REAL(NUMCRD),NUMCRD
421: 683     FORMAT(3(I6,1X),F8.4,1X,I5,' STUCT SNAP T T TID')420: 683     FORMAT(3(I6,1X),F8.4,1X,I5,' STUCT SNAP T T TID')
422: 421:         
423:         GLY_COUNT = 0422:         GLY_COUNT = 0
424: 423:         
425:         DO III = 1,NMRES424:         DO III = 1,NMRES
426:            IF (IRES(III).EQ.8) THEN425:            IF (IRES(III).EQ.8) THEN
427:               !!                pppcord(residue, xyz, numpro, atom types426:               !!                pppcord(residue, xyz, numpro, atom types
428:               PPPCORD(III, 1, 1, 1) = REAL(QMINP(J1,9*(III-1)+1- GLY_COUNT*3)) !  CA X427:               PPPCORD(III, 1, 1, 1) = REAL(QMINP(J1,9*(III-1)+1- GLY_COUNT*3)) !  CA X
429:               PPPCORD(III, 2, 1, 1) = REAL(QMINP(J1,9*(III-1)+2- GLY_COUNT*3)) !  CA Y428:               PPPCORD(III, 2, 1, 1) = REAL(QMINP(J1,9*(III-1)+2- GLY_COUNT*3)) !  CA Y
430:               PPPCORD(III, 3, 1, 1) = REAL(QMINP(J1,9*(III-1)+3- GLY_COUNT*3)) !  CA Z429:               PPPCORD(III, 3, 1, 1) = REAL(QMINP(J1,9*(III-1)+3- GLY_COUNT*3)) !  CA Z
431:               !    SWAP  CA for CB430:               !    SWAP  CA for CB
432:               PPPCORD(III, 1, 1, 2) = REAL(QMINP(j1,9*(III-1)+1- GLY_COUNT*3)) !  CB X431:               PPPCORD(III, 1, 1, 2) = REAL(QMINP(j1,9*(III-1)+1- GLY_COUNT*3)) !  CB X
433:               PPPCORD(III, 2, 1, 2) = REAL(QMINP(J1,9*(III-1)+2- GLY_COUNT*3)) !  CB Y432:               PPPCORD(III, 2, 1, 2) = REAL(QMINP(J1,9*(III-1)+2- GLY_COUNT*3)) !  CB Y
434:               PPPCORD(III, 3, 1, 2) = REAL(QMINP(J1,9*(III-1)+3- GLY_COUNT*3)) !  CB Z433:               PPPCORD(III, 3, 1, 2) = REAL(QMINP(J1,9*(III-1)+3- GLY_COUNT*3)) !  CB Z
435:               PPPCORD(III, 1, 1, 3) = REAL(QMINP(J1,9*(III-1)+4- GLY_COUNT*3)) !  O X434:               PPPCORD(III, 1, 1, 3) = REAL(QMINP(J1,9*(III-1)+4- GLY_COUNT*3)) !  O X
436:               PPPCORD(III, 2, 1, 3) = REAL(QMINP(J1,9*(III-1)+5- GLY_COUNT*3)) !  O Y435:               PPPCORD(III, 2, 1, 3) = REAL(QMINP(J1,9*(III-1)+5- GLY_COUNT*3)) !  O Y
437:               PPPCORD(III, 3, 1, 3) = REAL(QMINP(J1,9*(III-1)+6- GLY_COUNT*3)) !  O Z436:               PPPCORD(III, 3, 1, 3) = REAL(QMINP(J1,9*(III-1)+6- GLY_COUNT*3)) !  O Z
438: 437:               
439:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+1-GLY_COUNT*3),QMINP(J1,9*(III-1)+2-GLY_COUNT*3), &438:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+1-GLY_COUNT*3),QMINP(J1,9*(III-1)+2-GLY_COUNT*3), &
440:                    QMINP(J1,9*(III-1)+3-GLY_COUNT*3)439:                    QMINP(J1,9*(III-1)+3-GLY_COUNT*3)
441:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+1-GLY_COUNT*3),QMINP(J1,9*(III-1)+2-GLY_COUNT*3), &440:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+1-GLY_COUNT*3),QMINP(J1,9*(III-1)+2-GLY_COUNT*3), &
442:                    QMINP(J1,9*(III-1)+3-GLY_COUNT*3)441:                    QMINP(J1,9*(III-1)+3-GLY_COUNT*3)
443:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+4-GLY_COUNT*3),QMINP(J1,9*(III-1)+5-GLY_COUNT*3), &442:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+4-GLY_COUNT*3),QMINP(J1,9*(III-1)+5-GLY_COUNT*3), &
444:                    QMINP(J1,9*(III-1)+6-GLY_COUNT*3)443:                    QMINP(J1,9*(III-1)+6-GLY_COUNT*3)
445: 31            FORMAT('AM',3G25.15)444: 31            FORMAT('AM',3G25.15)
446: 445:               
447:               GLY_COUNT = GLY_COUNT +1446:               GLY_COUNT = GLY_COUNT +1
448:            ELSE447:            ELSE
449:               PPPCORD(III, 1, 1, 1) = REAL(QMINP(J1,9*(iii-1)+1- GLY_COUNT*3)) !  CA X448:               PPPCORD(III, 1, 1, 1) = REAL(QMINP(J1,9*(iii-1)+1- GLY_COUNT*3)) !  CA X
450:               PPPCORD(III, 2, 1, 1) = REAL(QMINP(J1,9*(III-1)+2- GLY_COUNT*3)) !  CA Y449:               PPPCORD(III, 2, 1, 1) = REAL(QMINP(J1,9*(III-1)+2- GLY_COUNT*3)) !  CA Y
451:               PPPCORD(III, 3, 1, 1) = REAL(QMINP(J1,9*(III-1)+3- GLY_COUNT*3)) !  CA Z450:               PPPCORD(III, 3, 1, 1) = REAL(QMINP(J1,9*(III-1)+3- GLY_COUNT*3)) !  CA Z
452:               PPPCORD(III, 1, 1, 2) = REAL(QMINP(J1,9*(III-1)+4- GLY_COUNT*3)) !  CB X451:               PPPCORD(III, 1, 1, 2) = REAL(QMINP(J1,9*(III-1)+4- GLY_COUNT*3)) !  CB X
453:               PPPCORD(III, 2, 1, 2) = REAL(QMINP(J1,9*(III-1)+5- GLY_COUNT*3)) !  CB Y452:               PPPCORD(III, 2, 1, 2) = REAL(QMINP(J1,9*(III-1)+5- GLY_COUNT*3)) !  CB Y
454:               PPPCORD(III, 3, 1, 2) = REAL(QMINP(J1,9*(III-1)+6- GLY_COUNT*3)) !  CB Z453:               PPPCORD(III, 3, 1, 2) = REAL(QMINP(J1,9*(III-1)+6- GLY_COUNT*3)) !  CB Z
455:               PPPCORD(III, 1, 1, 3) = REAL(QMINP(J1,9*(III-1)+7- GLY_COUNT*3)) !  O X454:               PPPCORD(III, 1, 1, 3) = REAL(QMINP(J1,9*(III-1)+7- GLY_COUNT*3)) !  O X
456:               PPPCORD(III, 2, 1, 3) = REAL(QMINP(J1,9*(III-1)+8- GLY_COUNT*3)) !  O Y455:               PPPCORD(III, 2, 1, 3) = REAL(QMINP(J1,9*(III-1)+8- GLY_COUNT*3)) !  O Y
457:               PPPCORD(III, 3, 1, 3) = REAL(QMINP(J1,9*(III-1)+9- GLY_COUNT*3)) !  O Z456:               PPPCORD(III, 3, 1, 3) = REAL(QMINP(J1,9*(III-1)+9- GLY_COUNT*3)) !  O Z
458: 457:               
459:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+1-GLY_COUNT*3),QMINP(J1,9*(III-1)+2-GLY_COUNT*3),&458:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+1-GLY_COUNT*3),QMINP(J1,9*(III-1)+2-GLY_COUNT*3),&
460:                    QMINP(J1,9*(III-1)+3-GLY_COUNT*3)459:                    QMINP(J1,9*(III-1)+3-GLY_COUNT*3)
461:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+4-GLY_COUNT*3),QMINP(J1,9*(III-1)+5-GLY_COUNT*3),&460:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+4-GLY_COUNT*3),QMINP(J1,9*(III-1)+5-GLY_COUNT*3),&
462:                    QMINP(J1,9*(III-1)+6-GLY_COUNT*3)461:                    QMINP(J1,9*(III-1)+6-GLY_COUNT*3)
463:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+7-GLY_COUNT*3),QMINP(J1,9*(III-1)+8-GLY_COUNT*3),&462:               WRITE(MYUNIT2,31) QMINP(J1,9*(III-1)+7-GLY_COUNT*3),QMINP(J1,9*(III-1)+8-GLY_COUNT*3),&
464:                    QMINP(J1,9*(III-1)+9-GLY_COUNT*3)463:                    QMINP(J1,9*(III-1)+9-GLY_COUNT*3)
465:            ENDIF464:            ENDIF
466:         ENDDO465:         ENDDO
467: 466:         
468:         DO III=1,NMRES467:         DO III=1,NMRES
469:            WRITE(AMHUNIT1,632)(PPPCORD(III,I3,1,1),I3=1,3),(PPPCORD(III,I3,1,2),I3=1,3),(PPPCORD(III,I3,1,3),I3=1,3)468:            WRITE(AMHUNIT1,632)(PPPCORD(III,I3,1,1),I3=1,3),(PPPCORD(III,I3,1,2),I3=1,3),(PPPCORD(III,I3,1,3),I3=1,3)
470: 632        FORMAT('CA: ',3(F8.3,1X),'CB: ',3(F8.3,1X),'OX: ', 3(F8.3,1X))469: 632        FORMAT('CA: ',3(F8.3,1X),'CB: ',3(F8.3,1X),'OX: ', 3(F8.3,1X))
471:            !632           FORMAT('CA: ',3(f25.15,1x),'CB: ',3(f25.15,1x),'Ox: ', 3(f25.15,1x))470:            !632           FORMAT('CA: ',3(f25.15,1x),'CB: ',3(f25.15,1x),'Ox: ', 3(f25.15,1x))
472:         END DO471:         END DO
473: 472:         
474:         DO III = 1+1, NMRES473:         DO III = 1+1, NMRES
475:            PPPCORD(III,1,1,4) = &474:            PPPCORD(III,1,1,4) = &
476:                 0.4831806D0*PPPCORD(III-1,1,1,1) + 0.7032820D0*PPPCORD(III,1,1,1) - 0.1864626D0*PPPCORD(III-1,1,1,3)475:                 0.4831806D0*PPPCORD(III-1,1,1,1) + 0.7032820D0*PPPCORD(III,1,1,1) - 0.1864626D0*PPPCORD(III-1,1,1,3)
477:            PPPCORD(III,2,1,4) = &476:            PPPCORD(III,2,1,4) = &
478:                 0.4831806D0*PPPCORD(III-1,2,1,1) + 0.7032820D0*PPPCORD(III,2,1,1) - 0.1864626D0*PPPCORD(III-1,2,1,3)477:                 0.4831806D0*PPPCORD(III-1,2,1,1) + 0.7032820D0*PPPCORD(III,2,1,1) - 0.1864626D0*PPPCORD(III-1,2,1,3)
479:            PPPCORD(III,3,1,4) = &478:            PPPCORD(III,3,1,4) = &
480:                 0.4831806D0*PPPCORD(III-1,3,1,1) + 0.7032820d0*PPPCORD(III,3,1,1) - 0.1864626D0*PPPCORD(III-1,3,1,3)479:                 0.4831806D0*PPPCORD(III-1,3,1,1) + 0.7032820d0*PPPCORD(III,3,1,1) - 0.1864626D0*PPPCORD(III-1,3,1,3)
481:         ENDDO480:         ENDDO
482: 481:         
483:         DO III = 1, NMRES-1482:         DO III = 1, NMRES-1
484:            PPPCORD(III,1,1,5) = &483:            PPPCORD(III,1,1,5) = &
485:                 0.4436538d0*PPPCORD(III,1,1,1)+0.2352006D0*PPPCORD(III+1,1,1,1)+0.3211455D0*PPPCORD(III,1,1,3)484:                 0.4436538d0*PPPCORD(III,1,1,1)+0.2352006D0*PPPCORD(III+1,1,1,1)+0.3211455D0*PPPCORD(III,1,1,3)
486:            PPPCORD(III,2,1,5) = &485:            PPPCORD(III,2,1,5) = &
487:                 0.4436538D0*PPPCORD(III,2,1,1)+0.2352006D0*PPPCORD(III+1,2,1,1)+0.3211455D0*PPPCORD(III,2,1,3)486:                 0.4436538D0*PPPCORD(III,2,1,1)+0.2352006D0*PPPCORD(III+1,2,1,1)+0.3211455D0*PPPCORD(III,2,1,3)
488:            PPPCORD(III,3,1,5) = &487:            PPPCORD(III,3,1,5) = &
489:                 0.4436538d0*PPPCORD(III,3,1,1)+0.2352006D0*PPPCORD(III+1,3,1,1)+0.3211455D0*PPPCORD(III,3,1,3)488:                 0.4436538d0*PPPCORD(III,3,1,1)+0.2352006D0*PPPCORD(III+1,3,1,1)+0.3211455D0*PPPCORD(III,3,1,3)
490:         ENDDO489:         ENDDO
491: 490:         
492:         DO III = 1, NMRES491:         DO III = 1, NMRES
493:            RES_TYPE = AMINOA(IRES(III))492:            RES_TYPE = AMINOA(IRES(III))
494: 493:            
495:            IF (III .NE. 1) THEN494:            IF (III .NE. 1) THEN
496:               ATOM_TYPE='N '495:               ATOM_TYPE='N '
497:               ID = 4496:               ID = 4
498:               WRITE(27,61)III,ATOM_TYPE,RES_TYPE,III,PPPCORD(III,1,1,ID),PPPCORD(III,2,1,ID),PPPCORD(III,3,1,ID),III497:               WRITE(27,61)III,ATOM_TYPE,RES_TYPE,III,PPPCORD(III,1,1,ID),PPPCORD(III,2,1,ID),PPPCORD(III,3,1,ID),III
499: 61            FORMAT('ATOM',4X,i3,2X,A2,2X,A3,3X,i3,4X,F8.3,F8.3,F8.3,2X,'1.00',2X,'0.00',6X,'TPDB',1x,I3)498: 61            FORMAT('ATOM',4X,i3,2X,A2,2X,A3,3X,i3,4X,F8.3,F8.3,F8.3,2X,'1.00',2X,'0.00',6X,'TPDB',1x,I3)
500: 499:               
501:            ENDIF500:            ENDIF
502: 501:            
503:            ATOM_TYPE='CA'502:            ATOM_TYPE='CA'
504:            ID = 1503:            ID = 1
505:            WRITE(27,61)III,ATOM_TYPE,RES_TYPE,III,PPPCORD(III,1,1,ID),PPPCORD(III,2,1,ID),PPPCORD(III,3,1,ID),III504:            WRITE(27,61)III,ATOM_TYPE,RES_TYPE,III,PPPCORD(III,1,1,ID),PPPCORD(III,2,1,ID),PPPCORD(III,3,1,ID),III
506: 505:            
507:            IF (RES_TYPE .NE. 'gly') THEN506:            IF (RES_TYPE .NE. 'gly') THEN
508:               ATOM_TYPE='CB'507:               ATOM_TYPE='CB'
509:               ID = 2508:               ID = 2
510:               WRITE(27,61)III,ATOM_TYPE,RES_TYPE,III,PPPCORD(III,1,1,ID),PPPCORD(III,2,1,ID),PPPCORD(III,3,1,ID),III509:               WRITE(27,61)III,ATOM_TYPE,RES_TYPE,III,PPPCORD(III,1,1,ID),PPPCORD(III,2,1,ID),PPPCORD(III,3,1,ID),III
511:            ENDIF510:            ENDIF
512: 511:            
513:            IF (III .NE. NMRES) THEN512:            IF (III .NE. NMRES) THEN
514:               ATOM_TYPE='C '513:               ATOM_TYPE='C '
515:               ID = 5514:               ID = 5
516:               WRITE(27,61)III,ATOM_TYPE,RES_TYPE,III,PPPCORD(III,1,1,ID),PPPCORD(III,2,1,ID),PPPCORD(III,3,1,ID),III515:               WRITE(27,61)III,ATOM_TYPE,RES_TYPE,III,PPPCORD(III,1,1,ID),PPPCORD(III,2,1,ID),PPPCORD(III,3,1,ID),III
517:            ENDIF516:            ENDIF
518: 517:            
519:            ATOM_TYPE='O '518:            ATOM_TYPE='O '
520:            ID = 3519:            ID = 3
521:            WRITE(27,61)III,ATOM_TYPE,RES_TYPE,III,PPPCORD(III,1,1,ID),PPPCORD(III,2,1,ID),PPPCORD(III,3,1,ID),III520:            WRITE(27,61)III,ATOM_TYPE,RES_TYPE,III,PPPCORD(III,1,1,ID),PPPCORD(III,2,1,ID),PPPCORD(III,3,1,ID),III
522: 521:            
523:         ENDDO522:         ENDDO
524:         CLOSE(27)523:         CLOSE(27)
525: 524:         
526:      ELSE IF (ARNO) THEN525:      ELSE IF (ARNO) THEN
527:         WRITE(MYUNIT2,'(A,F20.10)') 'N 0.0 0.0 ', 0.577D0526:         WRITE(MYUNIT2,'(A,F20.10)') 'N 0.0 0.0 ', 0.577D0
528:         WRITE(MYUNIT2,'(A,F20.10)') 'O 0.0 0.0 ',-0.577D0527:         WRITE(MYUNIT2,'(A,F20.10)') 'O 0.0 0.0 ',-0.577D0
529:         WRITE(MYUNIT2,65) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))528:         WRITE(MYUNIT2,65) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))
530: 65      FORMAT('AR ',3F20.10)529: 65      FORMAT('AR ',3F20.10)
531:      ELSE IF (TOSI.OR.WELCH) THEN530:      ELSE IF (TOSI.OR.WELCH) THEN
532:         DO J2=1,NATOMS531:         DO J2=1,NATOMS
533:            IF (ZSYM(J2).EQ.'PL') WRITE(MYUNIT2,'(A,3F20.10)') 'Na  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)532:            IF (ZSYM(J2).EQ.'PL') WRITE(MYUNIT2,'(A,3F20.10)') 'Na  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
534:            IF (ZSYM(J2).EQ.'MI') WRITE(MYUNIT2,'(A,3F20.10)') 'Cl  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)533:            IF (ZSYM(J2).EQ.'MI') WRITE(MYUNIT2,'(A,3F20.10)') 'Cl  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
535:         ENDDO534:         ENDDO
536:         ! hk286 - generalised Thomson problem, convert to Cartesians before printing535:         ! hk286 - generalised Thomson problem, convert to Cartesians before printing
537:      ELSE IF (GTHOMSONT) THEN536:      ELSE IF (GTHOMSONT) THEN
538:         CALL GTHOMSONANGTOC(DCOORDS(1:3*NATOMS), QMINP(J1, 1:3*NATOMS), NATOMS)537:         CALL GTHOMSONANGTOC(DCOORDS(1:3*NATOMS), QMINP(J1, 1:3*NATOMS), NATOMS) 
539:         IF (GTHOMPOT.EQ.7) THEN538:         IF (GTHOMPOT.EQ.7) THEN
540:            DO J2=1,NATOMS539:            DO J2=1,NATOMS
541:               IF (J2.LE.GTHOMSONBINN) THEN540:               IF (J2.LE.GTHOMSONBINN) THEN
542:                  WRITE(MYUNIT2,'(A,3F20.10)') 'LA  ',(DCOORDS(3*(J2-1)+J3),J3=1,3)541:                  WRITE(MYUNIT2,'(A,3F20.10)') 'LA  ',(DCOORDS(3*(J2-1)+J3),J3=1,3)
543:               ELSE542:               ELSE
544:                  WRITE(MYUNIT2,'(A,3F20.10)') 'LB  ',(DCOORDS(3*(J2-1)+J3),J3=1,3)543:                  WRITE(MYUNIT2,'(A,3F20.10)') 'LB  ',(DCOORDS(3*(J2-1)+J3),J3=1,3)
545:               ENDIF544:               ENDIF
546:            ENDDO545:            ENDDO
547:         ELSE546:         ELSE
548:            DO J2=1,NATOMS547:            DO J2=1,NATOMS
556:            ELSE555:            ELSE
557:               WRITE(MYUNIT2,'(A,3F20.10)') 'LB  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)556:               WRITE(MYUNIT2,'(A,3F20.10)') 'LB  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
558:            ENDIF557:            ENDIF
559:         ENDDO558:         ENDDO
560:      ELSE IF (GLJT) THEN559:      ELSE IF (GLJT) THEN
561:         J3=1560:         J3=1
562:         J4=1561:         J4=1
563:         DO J2=1,NATOMS562:         DO J2=1,NATOMS
564:            !563:            !
565:            DO J5=1,3564:            DO J5=1,3
566:               P(J5) = QMINP(J1,3*(J2-1)+J5)565:               P(J5) = QMINP(J1,3*(J2-1)+J5)                  
567:               IF(PERIODIC) THEN ! wrap back into the box566:               IF(PERIODIC) THEN ! wrap back into the box
568:                  P(J5)=P(J5)-BOX3D(J5)*ANINT(P(J5)/BOX3D(J5))567:                  P(J5)=P(J5)-BOX3D(J5)*ANINT(P(J5)/BOX3D(J5))
569:               ENDIF568:               ENDIF
570:            ENDDO569:            ENDDO
571:            !570:            !
572:            WRITE(ZSTR,'(I6)') J4571:            WRITE(ZSTR,'(I6)') J4
573:            ZSTR='L' // TRIM(ADJUSTL(ZSTR))572:            ZSTR='L' // TRIM(ADJUSTL(ZSTR))
574:            WRITE(MYUNIT2,'(A2,1X,3F20.10)') ZSTR,(P(J5),J5=1,3)573:            WRITE(MYUNIT2,'(A2,1X,3F20.10)') ZSTR,(P(J5),J5=1,3)
575:            J3=J3+1574:            J3=J3+1
576:            IF (J3.GT.NSPECIES(J4)) THEN575:            IF (J3.GT.NSPECIES(J4)) THEN
580:         ENDDO579:         ENDDO
581:         !580:         !
582:      ELSE IF (BGUPTAT .AND. .NOT. QALCST .AND. .NOT. MULTIPERMT) THEN581:      ELSE IF (BGUPTAT .AND. .NOT. QALCST .AND. .NOT. MULTIPERMT) THEN
583:         DO J2=1,NATOMS582:         DO J2=1,NATOMS
584:            IF (J2.LE.NTYPEA) THEN583:            IF (J2.LE.NTYPEA) THEN
585:               WRITE(MYUNIT2,'(A,3F20.10)') BGUPTANAME1, (QMINP(J1,3*(J2-1)+J3),J3=1,3)584:               WRITE(MYUNIT2,'(A,3F20.10)') BGUPTANAME1, (QMINP(J1,3*(J2-1)+J3),J3=1,3)
586:            ELSE585:            ELSE
587:               WRITE(MYUNIT2,'(A,3F20.10)') BGUPTANAME2, (QMINP(J1,3*(J2-1)+J3),J3=1,3)586:               WRITE(MYUNIT2,'(A,3F20.10)') BGUPTANAME2, (QMINP(J1,3*(J2-1)+J3),J3=1,3)
588:            ENDIF587:            ENDIF
589:         END DO588:         END DO
590: 589:         
591:      ELSE IF (AMBER.OR.MOLECULART) THEN590:      ELSE IF (AMBER.OR.MOLECULART) THEN
592:         DO J2=1,NATOMS591:         DO J2=1,NATOMS
593:            WRITE(MYUNIT2,'(A,3F20.10)') typech(J2)(1:1),(QMINP(J1,3*(J2-1)+J3),J3=1,3)592:            WRITE(MYUNIT2,'(A,3F20.10)') typech(J2)(1:1),(QMINP(J1,3*(J2-1)+J3),J3=1,3)
594:         ENDDO593:         ENDDO
595:      ELSE IF (AMBER12T) THEN594:      ELSE IF (AMBER12T) THEN
596:         !kr366> We use MYNODE + 1 throughout, need to use the same here!595:         !kr366> We use MYNODE + 1 throughout, need to use the same here!
597:         ! Create a string for J1596:         ! Create a string for J1
598:         WRITE(J1_STRING,'(I6)') J1597:         WRITE(J1_STRING,'(I6)') J1
599:         IF (DUMPSTRUCTURES) THEN598:         IF (DUMPSTRUCTURES) THEN
600:            IF (MPIT.OR.GBHT) THEN599:            IF (MPIT.OR.GBHT) THEN
623:                    &'.pdb', &622:                    &'.pdb', &
624:                    & LEN('coords.'//TRIM(ADJUSTL(J1_STRING))//'.pdb'))623:                    & LEN('coords.'//TRIM(ADJUSTL(J1_STRING))//'.pdb'))
625:            END IF624:            END IF
626:         END IF625:         END IF
627:         IF (MPIT.OR.GBHT) THEN626:         IF (MPIT.OR.GBHT) THEN
628:            WRITE(MYNODE_STRING,'(I6)') MYNODE + 1627:            WRITE(MYNODE_STRING,'(I6)') MYNODE + 1
629:            ! Pass the name, length of the string and whether or not to write a header.628:            ! Pass the name, length of the string and whether or not to write a header.
630:            CALL AMBER12_WRITE_XYZ(QMINP(J1,:), 'lowest.'//TRIM(ADJUSTL(MYNODE_STRING)), &629:            CALL AMBER12_WRITE_XYZ(QMINP(J1,:), 'lowest.'//TRIM(ADJUSTL(MYNODE_STRING)), &
631:                 LEN('lowest.'//TRIM(ADJUSTL(MYNODE_STRING))), &630:                 LEN('lowest.'//TRIM(ADJUSTL(MYNODE_STRING))), &
632:                 LOGICAL(.FALSE.,KIND=1))631:                 LOGICAL(.FALSE.,KIND=1))
633:         ELSE632:         ELSE 
634:            ! Pass the name, length of the string and whether or not to write a header.633:            ! Pass the name, length of the string and whether or not to write a header.
635:            CALL AMBER12_WRITE_XYZ(QMINP(J1,:), 'lowest', LEN('lowest'), LOGICAL(.FALSE.,KIND=1))634:            CALL AMBER12_WRITE_XYZ(QMINP(J1,:), 'lowest', LEN('lowest'), LOGICAL(.FALSE.,KIND=1))
636:         END IF635:         END IF
637:      ELSE IF (AMBERT) THEN636:      ELSE IF (AMBERT) THEN
638:         ! sf344> write out coordinates637:         ! sf344> write out coordinates
639:         COORDS1(1:3*NATOMS) = QMINP(J1,1:3*NATOMS)638:         COORDS1(1:3*NATOMS) = QMINP(J1,1:3*NATOMS)
640:         IF (DUMPSTRUCTURES) THEN639:         IF (DUMPSTRUCTURES) THEN
641:            WRITE(J1CHAR2,'(I3)') J1640:            WRITE(J1CHAR2,'(I3)') J1
642:            IF (MPIT.OR.GBHT) THEN641:            IF (MPIT.OR.GBHT) THEN
643:               CALL AMBERFINALIO(j1,MYUNIT2,MYNODE+1,tempstring,1,COORDS1(1:3*NATOMS))642:               CALL AMBERFINALIO(j1,MYUNIT2,MYNODE+1,tempstring,1,COORDS1(1:3*NATOMS))
653:            END IF652:            END IF
654:            DO J2=1,NATOMS653:            DO J2=1,NATOMS
655:               WRITE(MYUNIT3,'(3F28.20)') QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)654:               WRITE(MYUNIT3,'(3F28.20)') QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)
656:            ENDDO655:            ENDDO
657:            CLOSE(MYUNIT3)656:            CLOSE(MYUNIT3)
658:         ELSE657:         ELSE
659:            DO I1=1,NATOMS658:            DO I1=1,NATOMS
660:               WRITE(MYUNIT2,'(A2,3F20.10)') ih(m04+I1-1),COORDS1(3*I1-2),COORDS1(3*I1-1),COORDS1(3*I1)659:               WRITE(MYUNIT2,'(A2,3F20.10)') ih(m04+I1-1),COORDS1(3*I1-2),COORDS1(3*I1-1),COORDS1(3*I1)
661:            ENDDO660:            ENDDO
662:         ENDIF661:         ENDIF
663: 662:      
664:      ELSE IF (OPEPT) THEN663:      ELSE IF (OPEPT) THEN
665:         !we only have pdb files to write664:         !we only have pdb files to write
666:         WRITE(J1_STRING,'(I6)') J1665:         WRITE(J1_STRING,'(I6)') J1
667:         IF (MPIT.OR.GBHT) THEN666:         IF (MPIT.OR.GBHT) THEN
668:           WRITE(MYNODE_STRING,'(I6)') MYNODE + 1667:           WRITE(MYNODE_STRING,'(I6)') MYNODE + 1
669:           CALL OPEP_WRITE_PDB(NATOMS, QMINP(J1,:),'lowest.'//TRIM(ADJUSTL(J1_STRING))//&668:           CALL OPEP_WRITE_PDB(NATOMS, QMINP(J1,:),'lowest.'//TRIM(ADJUSTL(J1_STRING))//&
670:                    &'.'//TRIM(ADJUSTL(MYNODE_STRING))//'.pdb')669:                    &'.'//TRIM(ADJUSTL(MYNODE_STRING))//'.pdb')
671:           !also write start file to be used for OPTIM runs670:           !also write start file to be used for OPTIM runs
672:           MYUNIT3=GETUNIT()671:           MYUNIT3=GETUNIT()
673:           MYFILENAME2='start.'//TRIM(ADJUSTL(J1_STRING))//'.'//trim(adjustl(MYNODE_STRING))672:           MYFILENAME2='start.'//TRIM(ADJUSTL(J1_STRING))//'.'//trim(adjustl(MYNODE_STRING))
682:                     !also write start file to be used for OPTIM runs681:                     !also write start file to be used for OPTIM runs
683:           MYUNIT3=GETUNIT()682:           MYUNIT3=GETUNIT()
684:           MYFILENAME2='start.'//TRIM(ADJUSTL(J1_STRING))683:           MYFILENAME2='start.'//TRIM(ADJUSTL(J1_STRING))
685:           OPEN(MYUNIT3,FILE=trim(adjustl(MYFILENAME2)),STATUS="unknown",form="formatted")684:           OPEN(MYUNIT3,FILE=trim(adjustl(MYFILENAME2)),STATUS="unknown",form="formatted")
686:           DO J2=1,NATOMS685:           DO J2=1,NATOMS
687:               WRITE(MYUNIT3,'(3F28.20)') QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)686:               WRITE(MYUNIT3,'(3F28.20)') QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)
688:           ENDDO687:           ENDDO
689:           CLOSE(MYUNIT3)688:           CLOSE(MYUNIT3)
690: 689: 
691:         END IF690:         END IF
692: 691:         
693:      ELSE IF (CHIROT) THEN692:      ELSE IF (CHIROT) THEN
694:         CALL CHIRO_OUTPUT693:         CALL CHIRO_OUTPUT
695: 694:         
696:      ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT) THEN695:      ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT) THEN
697:         !696:         !
698:         ! determine the centre of coordinates, then centre them697:         ! determine the centre of coordinates, then centre them
699:         !698:         !
700:         SUMX=0.0D0699:         SUMX=0.0D0
701:         SUMY=0.0d0700:         SUMY=0.0d0
702:         SUMZ=0.0d0701:         SUMZ=0.0d0
703: 702:         
704:         DO J2=1,NATOMS/2703:         DO J2=1,NATOMS/2
705:            SUMX=SUMX+QMINP(J1,3*(J2-1)+1)704:            SUMX=SUMX+QMINP(J1,3*(J2-1)+1)
706:            SUMY=SUMY+QMINP(J1,3*(J2-1)+2)705:            SUMY=SUMY+QMINP(J1,3*(J2-1)+2)
707:            SUMZ=SUMZ+QMINP(J1,3*(J2-1)+3)706:            SUMZ=SUMZ+QMINP(J1,3*(J2-1)+3)
708:         ENDDO707:         ENDDO
709:         SUMX=2*SUMX/NATOMS708:         SUMX=2*SUMX/NATOMS
710:         SUMY=2*SUMY/NATOMS709:         SUMY=2*SUMY/NATOMS
711:         SUMZ=2*SUMZ/NATOMS710:         SUMZ=2*SUMZ/NATOMS
712:         DO J2=1,NATOMS/2711:         DO J2=1,NATOMS/2
713:            QMINP(J1,3*(J2-1)+1)=QMINP(J1,3*(J2-1)+1)-SUMX712:            QMINP(J1,3*(J2-1)+1)=QMINP(J1,3*(J2-1)+1)-SUMX
714:            QMINP(J1,3*(J2-1)+2)=QMINP(J1,3*(J2-1)+2)-SUMY713:            QMINP(J1,3*(J2-1)+2)=QMINP(J1,3*(J2-1)+2)-SUMY
715:            QMINP(J1,3*(J2-1)+3)=QMINP(J1,3*(J2-1)+3)-SUMZ714:            QMINP(J1,3*(J2-1)+3)=QMINP(J1,3*(J2-1)+3)-SUMZ
716:         ENDDO715:         ENDDO
717: 716:         
718:         DO j2=1,NATOMS/2717:         DO j2=1,NATOMS/2
719:            IF (PARAMONOVPBCX) THEN718:            IF (PARAMONOVPBCX) THEN
720:               ! ensure y component of particle 1 vector is within BoxLy/2 of zero.719:               ! ensure y component of particle 1 vector is within BoxLy/2 of zero. 
721:               ! If it isn't then subtract integer number of boxly's such that it is.720:               ! If it isn't then subtract integer number of boxly's such that it is.
722:               QMINP(J1,3*J2-2)=QMINP(J1,3*J2-2)-BOXLX*NINT(QMINP(J1,3*J2-2)/BOXLX)721:               QMINP(J1,3*J2-2)=QMINP(J1,3*J2-2)-BOXLX*NINT(QMINP(J1,3*J2-2)/BOXLX)
723:            ENDIF722:            ENDIF
724:            IF (PARAMONOVPBCY) THEN723:            IF (PARAMONOVPBCY) THEN
725:               ! ensure y component of particle 1 vector is within BoxLy/2 of zero.724:               ! ensure y component of particle 1 vector is within BoxLy/2 of zero. 
726:               ! If it isn't then subtract integer number of boxly's such that it is.725:               ! If it isn't then subtract integer number of boxly's such that it is.
727:               QMINP(J1,3*J2-1)=QMINP(J1,3*J2-1)-BOXLY*NINT(QMINP(J1,3*J2-1)/BOXLY)726:               QMINP(J1,3*J2-1)=QMINP(J1,3*J2-1)-BOXLY*NINT(QMINP(J1,3*J2-1)/BOXLY)
728:            ENDIF727:            ENDIF
729:            IF (PARAMONOVPBCZ) THEN728:            IF (PARAMONOVPBCZ) THEN
730:               ! ensure y component of particle 1 vector is within BoxLy/2 of zero.729:               ! ensure y component of particle 1 vector is within BoxLy/2 of zero. 
731:               ! If it isn't then subtract integer number of boxly's such that it is.730:               ! If it isn't then subtract integer number of boxly's such that it is.
732:               QMINP(J1,3*J2  )=QMINP(J1,3*J2  )-BOXLZ*NINT(QMINP(J1,3*J2  )/BOXLZ)731:               QMINP(J1,3*J2  )=QMINP(J1,3*J2  )-BOXLZ*NINT(QMINP(J1,3*J2  )/BOXLZ)
733:            ENDIF732:            ENDIF
734:         ENDDO733:         ENDDO
735: 734:         
736:         DO J2=1,NATOMS/2735:         DO J2=1,NATOMS/2
737:            WRITE(MYUNIT2,'(a5,2x,3f20.10,2x,a11,3f20.10)') 'H',QMINP(J1,3*(J2-1)+1), &736:            WRITE(MYUNIT2,'(a5,2x,3f20.10,2x,a11,3f20.10)') 'H',QMINP(J1,3*(J2-1)+1), &
738:                 QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3), &737:                 QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3), &
739:                 'atom_vector',QMINP(J1,3*NATOMS/2+3*(J2-1)+1), &738:                 'atom_vector',QMINP(J1,3*NATOMS/2+3*(J2-1)+1), &
740:                 QMINP(J1,3*NATOMS/2+3*(J2-1)+2),QMINP(J1,3*NATOMS/2+3*(J2-1)+3)739:                 QMINP(J1,3*NATOMS/2+3*(J2-1)+2),QMINP(J1,3*NATOMS/2+3*(J2-1)+3)
741:         ENDDO740:         ENDDO
742: 741:         
743:      ELSE IF (CHRMMT) THEN742:      ELSE IF (CHRMMT) THEN
744:         DO J2=1,NATOMS743:         DO J2=1,NATOMS
745:            WRITE(MYUNIT2,'(A,1X,3F20.10)') ZSYM(J2)(1:1),(QMINP(J1,3*(J2-1)+J3),J3=1,3)744:            WRITE(MYUNIT2,'(A,1X,3F20.10)') ZSYM(J2)(1:1),(QMINP(J1,3*(J2-1)+J3),J3=1,3)
746:         ENDDO745:         ENDDO
747:         !       csw34> This DO loop appeared to be be missing on 30/9/08 which would easily746:         !       csw34> This DO loop appeared to be be missing on 30/9/08 which would easily
748:         !       explain the output problems!747:         !       explain the output problems!
749:         DO J2=1,NATOMS748:         DO J2=1,NATOMS
750:            DCOORDS(3*(J2-1)+1)=QMINP(J1,3*(J2-1)+1)749:            DCOORDS(3*(J2-1)+1)=QMINP(J1,3*(J2-1)+1)
751:            DCOORDS(3*(J2-1)+2)=QMINP(J1,3*(J2-1)+2)750:            DCOORDS(3*(J2-1)+2)=QMINP(J1,3*(J2-1)+2)
752:            DCOORDS(3*(J2-1)+3)=QMINP(J1,3*(J2-1)+3)751:            DCOORDS(3*(J2-1)+3)=QMINP(J1,3*(J2-1)+3)
753:         ENDDO752:         ENDDO
754:         CALL CHARMMDUMP(DCOORDS,DBNAME(J1))753:         CALL CHARMMDUMP(DCOORDS,DBNAME(J1))
755: 754:         
756:         !    DC430 >755:         !    DC430 >
757:         !    |gd351> added patchy756:         !    |gd351> added patchy
758: 757:         
759:      ELSE IF (CAPBINT.OR.DBPT.OR.DBPTDT.OR.DMBLMT.OR.DMBLPYT.OR.LINRODT.OR.LWOTPT.OR.MSTBINT.OR.MSSTOCKT.OR.NCAPT.OR.NPAHT &758:      ELSE IF (CAPBINT.OR.DBPT.OR.DBPTDT.OR.DMBLMT.OR.DMBLPYT.OR.LINRODT.OR.LWOTPT.OR.MSTBINT.OR.MSSTOCKT.OR.NCAPT.OR.NPAHT &
760:           .OR. NTIPT .OR. STOCKAAT .OR. PAHAT .OR. PAHW99T .OR. TDHDT .OR. WATERDCT .OR. WATERKZT .OR. PATCHY .OR. PAPT   &759:           .OR. NTIPT .OR. STOCKAAT .OR. PAHAT .OR. PAHW99T .OR. TDHDT .OR. WATERDCT .OR. WATERKZT .OR. PATCHY .OR. PAPT   &
761:           .OR. PAPBINT .OR. PAPJANT .OR. POLYT .OR. PTSTSTT .OR. MORSEDPT) THEN760:           .OR. PAPBINT .OR. PAPJANT .OR. POLYT .OR. PTSTSTT .OR. MORSEDPT) THEN
762:         DO J2 = 1, NATOMS/2761:         DO J2 = 1, NATOMS/2
763:            WRITE(MYUNIT2,'(3f25.15)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)762:            WRITE(MYUNIT2,'(3f25.15)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)
764:         ENDDO763:         ENDDO
765:         DO J2 = 1, NATOMS/2764:         DO J2 = 1, NATOMS/2
766:            WRITE(MYUNIT2,'(3f25.15)') (QMINP(J1,3*NATOMS/2+3*(J2-1)+J3),J3=1,3)765:            WRITE(MYUNIT2,'(3f25.15)') (QMINP(J1,3*NATOMS/2+3*(J2-1)+J3),J3=1,3)
767:         ENDDO766:         ENDDO
768: 767:         
769:      ELSE IF (GBT.OR.GBDT.OR.GBDPT.OR.MSGBT) THEN768:      ELSE IF (GBT.OR.GBDT.OR.GBDPT.OR.MSGBT) THEN
770:         DO J2 = 1, NATOMS/2769:         DO J2 = 1, NATOMS/2
771:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)770:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)
772:         ENDDO771:         ENDDO
773:         DO J2 = 1, NATOMS/2772:         DO J2 = 1, NATOMS/2
774:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*NATOMS/2+3*(J2-1)+J3),J3=1,3)773:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*NATOMS/2+3*(J2-1)+J3),J3=1,3)
775:         ENDDO774:         ENDDO
776: 775:         
777:      ELSE IF (MLP3T.OR.MLPB3T.OR.MLQT.OR.MLPVB3T) THEN776:      ELSE IF (MLP3T.OR.MLPB3T.OR.MLQT.OR.MLPVB3T) THEN
778:         DO J2 = 1, NATOMS777:         DO J2 = 1, NATOMS
779:            WRITE(MYUNIT2,'(3G20.10)') QMINP(J1,J2)778:            WRITE(MYUNIT2,'(3G20.10)') QMINP(J1,J2)
780:         ENDDO779:         ENDDO
781: 780:         
782:      ELSE IF (GEMT) THEN781:      ELSE IF (GEMT) THEN
783:         DO J2 = 1, NATOMS782:         DO J2 = 1, NATOMS
784:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)783:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)
785:         ENDDO784:         ENDDO
786: 785:         
787:      ELSE IF (BLNT.AND.(.NOT.P46).AND.(.NOT.G46)) THEN786:      ELSE IF (BLNT.AND.(.NOT.P46).AND.(.NOT.G46)) THEN
788:         !787:         !
789:         ! this writes 'lowest' in xyz (Xmakemol) format788:         ! this writes 'lowest' in xyz (Xmakemol) format
790:         !789:         !
791:         WRITE(MYUNIT,'(A,I6,A)') ' in finalio BLN block MYUNIT2=',MYUNIT2790:         WRITE(MYUNIT,'(A,I6,A)') ' in finalio BLN block MYUNIT2=',MYUNIT2
792:         DO J2=1,NATOMS791:         DO J2=1,NATOMS
793:            WRITE(MYUNIT2,'(2A1,1X,3F20.10)') BEADLETTER(J2),'L',(QMINP(J1,3*(J2-1)+J3),J3=1,3)792:            WRITE(MYUNIT2,'(2A1,1X,3F20.10)') BEADLETTER(J2),'L',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
794:         ENDDO793:         ENDDO
795:      ELSE IF(QALCST.OR.MULTIPERMT.OR.MLJT.OR.MGUPTAT.OR.MSCT) THEN794:      ELSE IF(QALCST.OR.MULTIPERMT.OR.MLJT.OR.MGUPTAT.OR.MSCT) THEN
796:         !795:         !
797:         CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)796:         CALL SET_ATOMLISTS(QMINT(J1,1:NATOMS),1)
798:         CALL SET_LABELS(QMINT(J1,1:NATOMS),1)797:         CALL SET_LABELS(QMINT(J1,1:NATOMS),1)
799:         IF(SPECMASST) THEN798:         IF(SPECMASST) THEN 
800:            EDUMMY = SUM(ATMASS)799:            EDUMMY = SUM(ATMASS)
801:            WRITE(MYUNIT,'(A,I6,A,E20.10)') &800:            WRITE(MYUNIT,'(A,I6,A,E20.10)') &
802:                 'finalio> Min ',J1,' has total mass ', EDUMMY801:                 'finalio> Min ',J1,' has total mass ', EDUMMY
803:         ENDIF802:         ENDIF
804:         !803:         !
805:         IF(STRESST) CALL CALC_STRESS(QMINP(J1,1:3*NATOMS),EDUMMY)804:         IF(STRESST) CALL CALC_STRESS(QMINP(J1,1:3*NATOMS),EDUMMY)
806:         !805:         !
807:         DO J2=1,NSPECIES(0) !ds656> Should not exceed 10806:         DO J2=1,NSPECIES(0) !ds656> Should not exceed 10
808:            IF(SPECLABELST) THEN807:            IF(SPECLABELST) THEN
809:               ATOM_TYPE = SPECLABELS(J2)808:               ATOM_TYPE = SPECLABELS(J2)
810:            ELSE809:            ELSE
811:               WRITE(ATOM_TYPE,'(I1)') J2810:               WRITE(ATOM_TYPE,'(I1)') J2
812:               ATOM_TYPE='L' // TRIM(ADJUSTL(ATOM_TYPE))811:               ATOM_TYPE='L' // TRIM(ADJUSTL(ATOM_TYPE))
813:            ENDIF812:            ENDIF
814:            DO J3=1,ATOMLISTS(J2,1,0)813:            DO J3=1,ATOMLISTS(J2,1,0)
815:               J4=ATOMLISTS(J2,1,J3) ! Actual atom index814:               J4=ATOMLISTS(J2,1,J3) ! Actual atom index
816:               DO J5=1,3815:               DO J5=1,3
817:                  P(J5) = QMINP(J1,3*(J4-1)+J5)816:                  P(J5) = QMINP(J1,3*(J4-1)+J5)                  
818:                  IF(PERIODIC) THEN ! wrap back into the box817:                  IF(PERIODIC) THEN ! wrap back into the box
819:                     P(J5)=P(J5)-BOX3D(J5)*ANINT(P(J5)/BOX3D(J5))818:                     P(J5)=P(J5)-BOX3D(J5)*ANINT(P(J5)/BOX3D(J5))
820:                  ENDIF819:                  ENDIF
821:               ENDDO820:               ENDDO
822:               WRITE(MYUNIT2,'(A,3(1X,F20.10))',ADVANCE='NO') &821:               WRITE(MYUNIT2,'(A,3(1X,F20.10))',ADVANCE='NO') &
823:                    ATOM_TYPE,(P(J5), J5=1,3)822:                    ATOM_TYPE,(P(J5), J5=1,3)
824:               IF(STRESST) THEN ! Tack on local stresses823:               IF(STRESST) THEN ! Tack on local stresses
825:                  IF(STRESS_MODE==2) THEN824:                  IF(STRESS_MODE==2) THEN
826:                     ! Just print the uniqe elements of the stress825:                     ! Just print the uniqe elements of the stress
827:                     ! tensor826:                     ! tensor
848:      ELSE847:      ELSE
849:         IF (CSMT.AND.(.NOT.SYMMETRIZECSM)) THEN848:         IF (CSMT.AND.(.NOT.SYMMETRIZECSM)) THEN
850:            WRITE(MYUNIT3,'(I6)') NATOMS849:            WRITE(MYUNIT3,'(I6)') NATOMS
851:            WRITE(MYUNIT3,'(A,I6,2(A,G20.10))') 'averaged structure for final solution ',J1, &850:            WRITE(MYUNIT3,'(A,I6,2(A,G20.10))') 'averaged structure for final solution ',J1, &
852:                 ' CSM=',QMINAV(J1),' CSM for reference structure=',QMIN(J1)851:                 ' CSM=',QMINAV(J1),' CSM for reference structure=',QMIN(J1)
853:            WRITE(MYUNIT3,30) (QMINPCSMAV(J1,J2),J2=1,3*(NATOMS-NS))852:            WRITE(MYUNIT3,30) (QMINPCSMAV(J1,J2),J2=1,3*(NATOMS-NS))
854:         ENDIF853:         ENDIF
855:         WRITE(MYUNIT2,30) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))854:         WRITE(MYUNIT2,30) (QMINP(J1,J2),J2=1,3*(NATOMS-NS))
856: 30      FORMAT('LA ',3F20.10)855: 30      FORMAT('LA ',3F20.10)
857:      ENDIF856:      ENDIF
858: 857:      
859:      !|gd351>858:      !|gd351>
860:      IF (ASAOOS) THEN859:      IF (ASAOOS) THEN
861:         LUNIT=GETUNIT()860:         LUNIT=GETUNIT()
862:         OPEN(LUNIT,file='particles.xyz')861:         OPEN(LUNIT,file='particles.xyz')
863:         WRITE(LUNIT,*) NATOMS862:         WRITE(LUNIT,*) NATOMS
864:         WRITE(LUNIT,*) ' '863:         WRITE(LUNIT,*) ' '
865:         DO J2=1,NATOMS864:         DO J2=1,NATOMS
866:            WRITE(LUNIT,'(A,3F20.10)') 'H  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)865:            WRITE(LUNIT,'(A,3F20.10)') 'H  ',(QMINP(J1,3*(J2-1)+J3),J3=1,3)
867:         ENDDO866:         ENDDO
868:         CLOSE(LUNIT)867:         CLOSE(LUNIT)
869:      END IF868:      END IF
870:      !<gd351|869:      !<gd351|
871: 870:      
872:      IF ((NS.GT.0).AND.(.NOT.(WELCH.OR.TOSI))) THEN871:      IF ((NS.GT.0).AND.(.NOT.(WELCH.OR.TOSI))) THEN
873:         IF (MSORIGT.OR.FRAUSIT) THEN872:         IF (MSORIGT.OR.FRAUSIT) THEN
874:            WRITE(MYUNIT2,40) (QMINP(J1,J2),J2=3*(NATOMS-NS)+1,3*NATOMS)873:            WRITE(MYUNIT2,40) (QMINP(J1,J2),J2=3*(NATOMS-NS)+1,3*NATOMS)
875: 40         FORMAT('Si',3F20.10)874: 40         FORMAT('Si',3F20.10)
876:         ELSE IF (MSTRANST) THEN875:         ELSE IF (MSTRANST) THEN
877:            WRITE(MYUNIT2,40) (QMINP(J1,J2),J2=3*(NATOMS-NS)+1,3*NATOMS)876:            WRITE(MYUNIT2,40) (QMINP(J1,J2),J2=3*(NATOMS-NS)+1,3*NATOMS)
878:         ELSE877:         ELSE
879:            WRITE(MYUNIT2,50) (QMINP(J1,J2),J2=3*(NATOMS-NS)+1,3*NATOMS)878:            WRITE(MYUNIT2,50) (QMINP(J1,J2),J2=3*(NATOMS-NS)+1,3*NATOMS)
880: 50         FORMAT('LB',3F20.10)879: 50         FORMAT('LB',3F20.10)
881:         ENDIF880:         ENDIF
883:      IF (AMBER) CALL AMBERDUMP(J1,QMINP)882:      IF (AMBER) CALL AMBERDUMP(J1,QMINP)
884: !883: !
885: ! Reorder the addressable particles into clusters.884: ! Reorder the addressable particles into clusters.
886: !885: !
887:      IF (REORDERADDT) THEN886:      IF (REORDERADDT) THEN
888:         WRITE(REORDERUNIT,'(I10)') NATOMS887:         WRITE(REORDERUNIT,'(I10)') NATOMS
889:         WRITE(REORDERUNIT,10) J1, QMIN(J1), FF(J1), NPCALL_QMIN(J1)888:         WRITE(REORDERUNIT,10) J1, QMIN(J1), FF(J1), NPCALL_QMIN(J1)
890:         ADDDONE(1:NATOMS)=.FALSE.889:         ADDDONE(1:NATOMS)=.FALSE.
891:         DO J2=1,NATOMS/NADDTARGET890:         DO J2=1,NATOMS/NADDTARGET
892:            MJ2=(J2-1)*NADDTARGET+1891:            MJ2=(J2-1)*NADDTARGET+1
893: !           CLOSEST(1)=MJ2892:            CLOSEST(1)=MJ2
894:            WRITE(REORDERUNIT,'(A2,1X,3G20.10)') 'LA ',QMINP(J1,3*(MJ2-1)+1),QMINP(J1,3*(MJ2-1)+2),QMINP(J1,3*(MJ2-1)+3)893:            WRITE(REORDERUNIT,'(A2,1X,3G20.10)') 'LA ',QMINP(J1,3*(MJ2-1)+1),QMINP(J1,3*(MJ2-1)+2),QMINP(J1,3*(MJ2-1)+3)
895:            ADDDONE(MJ2)=.TRUE.894:            ADDDONE(MJ2)=.TRUE.
896:            DO J3=2,NADDTARGET895:            DO J3=2,NADDTARGET
897:               DMIN=1.0D100896:               DMIN=1.0D100
898:               DO J4=1,NATOMS/NADDTARGET897:               DO J4=1,NATOMS/NADDTARGET
899:                  IF (ADDDONE((J4-1)*NADDTARGET+J3)) CYCLE898:                  IF (ADDDONE((J4-1)*NADDTARGET+J3)) CYCLE
900:                  DIST=(QMINP(J1,3*(MJ2-1)+1)-QMINP(J1,3*((J4-1)*NADDTARGET+J3-1)+1))**2+ &899:                  DIST=(QMINP(J1,3*(MJ2-1)+1)-QMINP(J1,3*((J4-1)*NADDTARGET+J3-1)+1))**2+ &
901:   &                   (QMINP(J1,3*(MJ2-1)+2)-QMINP(J1,3*((J4-1)*NADDTARGET+J3-1)+2))**2+ &900:   &                   (QMINP(J1,3*(MJ2-1)+2)-QMINP(J1,3*((J4-1)*NADDTARGET+J3-1)+2))**2+ &
902:   &                   (QMINP(J1,3*(MJ2-1)+3)-QMINP(J1,3*((J4-1)*NADDTARGET+J3-1)+3))**2901:   &                   (QMINP(J1,3*(MJ2-1)+3)-QMINP(J1,3*((J4-1)*NADDTARGET+J3-1)+3))**2
903:                  IF (DIST.LT.DMIN) THEN902:                  IF (DIST.LT.DMIN) THEN
904:                     DMIN=DIST903:                     DMIN=DIST
905:                     NCLOSE=(J4-1)*NADDTARGET+J3904:                     NCLOSE=(J4-1)*NADDTARGET+J3
906:                  ENDIF905:                  ENDIF
907:                  WRITE(MYUNIT,'(A,4I6,2G20.10)') 'J2,MJ2,J3,NCLOSE,DIST,DMIN=',J2,MJ2,J3,NCLOSE,DIST,DMIN906:                  WRITE(MYUNIT,'(A,4I6,2G20.10)') 'J2,MJ2,J3,NCLOSE,DIST,DMIN=',J2,MJ2,J3,NCLOSE,DIST,DMIN
908:               ENDDO907:               ENDDO
909: !              CLOSEST(J3)=NCLOSE908:               CLOSEST(J3)=NCLOSE
910:               ADDDONE(NCLOSE)=.TRUE.909:               ADDDONE(NCLOSE)=.TRUE.
911:               WRITE(MYUNIT,'(A,4I6)') 'J2,MJ2,J3, final closest=',J2,MJ2,J3,NCLOSE910:               WRITE(MYUNIT,'(A,4I6)') 'J2,MJ2,J3, final closest=',J2,MJ2,J3,NCLOSE
912:               WRITE(REORDERUNIT,'(A2,1X,3G20.10)') 'LA ',QMINP(J1,3*(NCLOSE-1)+1),QMINP(J1,3*(NCLOSE-1)+2),QMINP(J1,3*(NCLOSE-1)+3)911:               WRITE(REORDERUNIT,'(A2,1X,3G20.10)') 'LA ',QMINP(J1,3*(NCLOSE-1)+1),QMINP(J1,3*(NCLOSE-1)+2),QMINP(J1,3*(NCLOSE-1)+3)
913:            ENDDO912:            ENDDO
914:         ENDDO913:         ENDDO
915:      ENDIF914:      ENDIF
916:      IF (LJADD3T.AND.(NADDTARGET.LE.26)) THEN915:      IF (LJADD3T.AND.(NADDTARGET.LE.26)) THEN
917:         WRITE(LJGOUNIT,'(I10)') NATOMS916:         WRITE(LJGOUNIT,'(I10)') NATOMS
918:         WRITE(LJGOUNIT,10) J1, QMIN(J1), FF(J1), NPCALL_QMIN(J1)917:         WRITE(LJGOUNIT,10) J1, QMIN(J1), FF(J1), NPCALL_QMIN(J1)
919:         DO J2=1,NATOMS918:         DO J2=1,NATOMS
920:            MJ2=MOD(J2-1,NADDTARGET)+1919:            MJ2=MOD(J2-1,NADDTARGET)+1
921:            WRITE(LJGOUNIT,'(A2,1X,3G20.10)') ALPHABET(MJ2),QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)920:            WRITE(LJGOUNIT,'(A2,1X,3G20.10)') ALPHABET(MJ2),QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)
922:         ENDDO921:         ENDDO
923:      ENDIF922:      ENDIF
924:      !923:      !
925:   ENDDO savemin924:   ENDDO savemin
926:   !925:   !
927:   ! End of loop over dump to file lowest or equivalent.926:   ! End of loop over dump to file lowest or equivalent. 
928:   !927:   !
929:   CLOSE(MYUNIT2)928:   CLOSE(MYUNIT2)
930:   IF (CSMT.AND.(.NOT.SYMMETRIZECSM)) CLOSE(MYUNIT3)929:   IF (CSMT.AND.(.NOT.SYMMETRIZECSM)) CLOSE(MYUNIT3)
931:   !930:   !
932:   !     csw34> New loop for dumping interaction energy files if A9INTE is specified931:   !     csw34> New loop for dumping interaction energy files if A9INTE is specified
933:   !     Added the missing IF block to test for A9INTE 9/12/09 DJW932:   !     Added the missing IF block to test for A9INTE 9/12/09 DJW
934:   !933:   !
935:   IF (A9INTET) THEN934:   IF (A9INTET) THEN
936:      IF (MPIT.OR.GBHT) THEN935:      IF (MPIT.OR.GBHT) THEN
937:         WRITE (ISTR, '(I10)') MYUNIT-22980+1936:         WRITE (ISTR, '(I10)') MYUNIT-22980+1
951:         WRITE(MYUNIT2,*) NATOMS950:         WRITE(MYUNIT2,*) NATOMS
952:         !     csw34> write header to intelowest for current minimum951:         !     csw34> write header to intelowest for current minimum
953:         WRITE(MYUNIT2,10) J1, INTEQMIN(J1), INTEFF(J1)952:         WRITE(MYUNIT2,10) J1, INTEQMIN(J1), INTEFF(J1)
954:         !     sf344> write out coordinates953:         !     sf344> write out coordinates
955:         COORDS1(1:3*NATOMS) = INTEQMINP(J1,1:3*NATOMS)954:         COORDS1(1:3*NATOMS) = INTEQMINP(J1,1:3*NATOMS)
956:         IF (DUMPSTRUCTURES) THEN955:         IF (DUMPSTRUCTURES) THEN
957:            CALL INTEFINALIO(j1,MYUNIT2,AMBFINALIO_NODE,'0',0,COORDS1(1:3*NATOMS))956:            CALL INTEFINALIO(j1,MYUNIT2,AMBFINALIO_NODE,'0',0,COORDS1(1:3*NATOMS))
958:            WRITE(J1CHAR2,'(I3)') J1957:            WRITE(J1CHAR2,'(I3)') J1
959:            WRITE(J1CHAR,'(A,A)') 'intecoords.',TRIM(ADJUSTL(J1CHAR2))958:            WRITE(J1CHAR,'(A,A)') 'intecoords.',TRIM(ADJUSTL(J1CHAR2))
960:            OPEN(UNIT=226,FILE=trim(adjustl(J1CHAR)),STATUS='UNKNOWN')959:            OPEN(UNIT=226,FILE=trim(adjustl(J1CHAR)),STATUS='UNKNOWN')
961: 960:            
962:            DO J2=1,NATOMS961:            DO J2=1,NATOMS
963:               WRITE(226,'(3F28.20)') INTEQMINP(J1,3*(J2-1)+1),INTEQMINP(J1,3*(J2-1)+2),INTEQMINP(J1,3*(J2-1)+3)962:               WRITE(226,'(3F28.20)') INTEQMINP(J1,3*(J2-1)+1),INTEQMINP(J1,3*(J2-1)+2),INTEQMINP(J1,3*(J2-1)+3)
964:            ENDDO963:            ENDDO
965:            CLOSE(226)964:            CLOSE(226)
966:            WRITE(J1CHAR2,'(I3)') J1965:            WRITE(J1CHAR2,'(I3)') J1
967:            WRITE(J1CHAR,'(A,A)') 'intecoords.',TRIM(ADJUSTL(J1CHAR2))966:            WRITE(J1CHAR,'(A,A)') 'intecoords.',TRIM(ADJUSTL(J1CHAR2))
968:            OPEN(UNIT=226,FILE=trim(adjustl(J1CHAR)),STATUS='UNKNOWN')967:            OPEN(UNIT=226,FILE=trim(adjustl(J1CHAR)),STATUS='UNKNOWN')
969: 968:            
970:            DO J2=1,NATOMS969:            DO J2=1,NATOMS
971:               WRITE(226,'(3F28.20)') INTEQMINP(J1,3*(J2-1)+1),INTEQMINP(J1,3*(J2-1)+2),INTEQMINP(J1,3*(J2-1)+3)970:               WRITE(226,'(3F28.20)') INTEQMINP(J1,3*(J2-1)+1),INTEQMINP(J1,3*(J2-1)+2),INTEQMINP(J1,3*(J2-1)+3)
972:            ENDDO971:            ENDDO
973:            CLOSE(226)972:            CLOSE(226)
974: 973:            
975:         ELSE974:         ELSE
976:            DO I1=1,NATOMS975:            DO I1=1,NATOMS
977:               WRITE(MYUNIT2,'(A2,3F20.10)') ih(m04+I1-1),COORDS1(3*I1-2),COORDS1(3*I1-1),COORDS1(3*I1)976:               WRITE(MYUNIT2,'(A2,3F20.10)') ih(m04+I1-1),COORDS1(3*I1-2),COORDS1(3*I1-1),COORDS1(3*I1)
978:            ENDDO977:            ENDDO
979:         ENDIF978:         ENDIF
980:      ENDDO979:      ENDDO
981:   ENDIF980:   ENDIF
982:   !981:   !
983:   !  End of loop over dump to file intelowest982:   !  End of loop over dump to file intelowest
984:   !983:   !        
985:   ! csw34> Output for the HBONDMATRIX keyword (subroutine in hbondmatrix.f90)984:   ! csw34> Output for the HBONDMATRIX keyword (subroutine in hbondmatrix.f90)
986:   IF ((AMBERT.OR.AMBER12T.OR.(CUDAT.AND.(CUDAPOT.EQ.'A'))).AND.HBONDMATRIX) CALL HBONDMATRIXFINALIO()985:   IF ((AMBERT.OR.AMBER12T.OR.(CUDAT.AND.(CUDAPOT.EQ.'A'))).AND.HBONDMATRIX) CALL HBONDMATRIXFINALIO()
987: 986:   
988:   CLOSE(MYUNIT2)987:   CLOSE(MYUNIT2)
989: 988:   
990:   !     csw34> Edits to the RMS keyword989:   !     csw34> Edits to the RMS keyword
991:   IF (CHRMMT.AND.RMST) THEN990:   IF (CHRMMT.AND.RMST) THEN
992:      !        IF (PROGRESS) THEN991:      !        IF (PROGRESS) THEN
993:      !           DCOORDS(1:3*NATOMS)=RMSCOOR(1,1:3*NATOMS)992:      !           DCOORDS(1:3*NATOMS)=RMSCOOR(1,1:3*NATOMS)
994:      !           IF(RMSBEST(1,2)<0.D0) CALL CHARMMDUMP(DCOORDS,'closestrms')993:      !           IF(RMSBEST(1,2)<0.D0) CALL CHARMMDUMP(DCOORDS,'closestrms')
995:      !           WRITE(MYUNIT,'(A9,F8.5)') 'RMSDmin= ',RMSBEST(1,1)994:      !           WRITE(MYUNIT,'(A9,F8.5)') 'RMSDmin= ',RMSBEST(1,1)
996:      !        ELSE995:      !        ELSE
997:      ! REMEMBER TO RE-INDENT THE BELOW IF UNCOMMENTING ABOVE!996:      ! REMEMBER TO RE-INDENT THE BELOW IF UNCOMMENTING ABOVE!
998:      OPEN(UNIT=MYUNIT2,FILE='rmsbest.'//TRIM(ADJUSTL(ISTR)),STATUS='UNKNOWN')997:      OPEN(UNIT=MYUNIT2,FILE='rmsbest.'//TRIM(ADJUSTL(ISTR)),STATUS='UNKNOWN')
999:      DO J2=1,RMSSAVE998:      DO J2=1,RMSSAVE
1000:         WRITE(MYUNIT2,'(I6,F6.3,F15.5)')J2,RMSBEST(J2,1),RMSBEST(J2,2)999:         WRITE(MYUNIT2,'(I6,F6.3,F15.5)')J2,RMSBEST(J2,1),RMSBEST(J2,2)
1001:         WRITE(CRMS,'(I6)') J21000:         WRITE(CRMS,'(I6)') J2
1002:         DCOORDS(1:3*NATOMS)=RMSCOOR(J2,1:3*NATOMS)1001:         DCOORDS(1:3*NATOMS)=RMSCOOR(J2,1:3*NATOMS)
1003:         IF(RMSBEST(J2,2)<0.D0) CALL CHARMMDUMP(DCOORDS,'rms.'//TRIM(ADJUSTL(ISTR))//'.'//TRIM(ADJUSTL(CRMS)))1002:         IF(RMSBEST(J2,2)<0.D0) CALL CHARMMDUMP(DCOORDS,'rms.'//TRIM(ADJUSTL(ISTR))//'.'//TRIM(ADJUSTL(CRMS)))
1004:      ENDDO1003:      ENDDO
1005:      CLOSE(MYUNIT2)1004:      CLOSE(MYUNIT2)
1006:      !        ENDIF1005:      !        ENDIF
1007:   ENDIF1006:   ENDIF
1008: 1007:   
1009:   IF (LJCOULT) THEN1008:   IF (LJCOULT) THEN
1010:      OPEN(UNIT=AMHUNIT1,FILE='ljcoul.xyz',STATUS='UNKNOWN')1009:      OPEN(UNIT=AMHUNIT1,FILE='ljcoul.xyz',STATUS='UNKNOWN')
1011:      DO J1=1,NSAVE1010:      DO J1=1,NSAVE
1012:         WRITE(AMHUNIT1,'(I6)') NATOMS1011:         WRITE(AMHUNIT1,'(I6)') NATOMS
1013:         WRITE(AMHUNIT1,10) J1, QMIN(J1), FF(J1)1012:         WRITE(AMHUNIT1,10) J1, QMIN(J1), FF(J1)
1014:         DO J2=1,NATOMS1013:         DO J2=1,NATOMS
1015:            !              Use "O" atom type to highlight charged particles and "N" for the neutral ones.1014:            !              Use "O" atom type to highlight charged particles and "N" for the neutral ones.
1016:            IF (J2.LE.COULN) THEN1015:            IF (J2.LE.COULN) THEN
1017:               WRITE(AMHUNIT1,'(A4,3F18.10,A12,3F18.10)') 'O ',&1016:               WRITE(AMHUNIT1,'(A4,3F18.10,A12,3F18.10)') 'O ',&
1018:                    QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)1017:                    QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)
1019:            ELSE1018:            ELSE
1020:               WRITE(AMHUNIT1,'(A4,3F18.10,A12,3F18.10)') 'N ',&1019:               WRITE(AMHUNIT1,'(A4,3F18.10,A12,3F18.10)') 'N ',&
1021:                    QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)1020:                    QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)
1022:            END IF1021:            END IF
1023:         ENDDO1022:         ENDDO
1024:      ENDDO1023:      ENDDO
1025:      CLOSE(AMHUNIT1)1024:      CLOSE(AMHUNIT1)
1026: 1025:      
1027: 1026:      
1028:   ELSE IF (STOCKT) THEN1027:   ELSE IF (STOCKT) THEN
1029:      LUNIT=GETUNIT()1028:      LUNIT=GETUNIT()
1030:      OPEN(UNIT=LUNIT,FILE='stock.xyz',STATUS='UNKNOWN')1029:      OPEN(UNIT=LUNIT,FILE='stock.xyz',STATUS='UNKNOWN')
1031:      DO J1=1,NSAVE1030:      DO J1=1,NSAVE
1032:         WRITE(LUNIT,'(I6)') NATOMS/21031:         WRITE(LUNIT,'(I6)') NATOMS/2
1033:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1032:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1034:         DO J2=1,NATOMS/21033:         DO J2=1,NATOMS/2
1035:            WRITE (LUNIT,'(A4,3F18.10,A12,3F18.10)') 'LA ', &1034:            WRITE (LUNIT,'(A4,3F18.10,A12,3F18.10)') 'LA ', &
1036:                 QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3), &1035:                 QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3), &
1037:                 'atom_vector', &1036:                 'atom_vector', &
1038:                 SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))*COS(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+2)), &1037:                 SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))*COS(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+2)), &
1039:                 SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))*SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+2)), &1038:                 SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))*SIN(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+2)), &
1040:                 COS(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))1039:                 COS(QMINP(J1,3*(NATOMS/2)+3*(J2-1)+1))
1041:         ENDDO1040:         ENDDO
1042:      ENDDO1041:      ENDDO
1043:      CLOSE(LUNIT)1042:      CLOSE(LUNIT)
1044: 1043:      
1045:   ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT.OR.PYT) THEN1044:   ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT.OR.PYT) THEN
1046:      DO J1=1,NSAVE1045:      DO J1=1,NSAVE
1047:         WRITE(J1CHAR2,'(I3)') J11046:         WRITE(J1CHAR2,'(I3)') J1
1048:         IF (MPIT.OR.GBHT) THEN1047:         IF (MPIT.OR.GBHT) THEN
1049:            WRITE (ISTR, '(I10)') MYNODE+11048:            WRITE (ISTR, '(I10)') MYNODE+1
1050:            MYUNIT2=GETUNIT()1049:            MYUNIT2=GETUNIT()
1051:            MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))//'.'//trim(adjustl(istr))1050:            MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))//'.'//trim(adjustl(istr))
1052:            OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")1051:            OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")
1053:         ELSE1052:         ELSE
1054:            MYUNIT2=GETUNIT()1053:            MYUNIT2=GETUNIT()
1055:            MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))1054:            MYFILENAME2='coords.'//TRIM(ADJUSTL(J1CHAR2))
1056:            OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")1055:            OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")
1057:         END IF1056:         END IF
1058: 1057:         
1059:         DO J2=1,NATOMS1058:         DO J2=1,NATOMS
1060:            WRITE(MYUNIT2,'(3F28.20)') QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)1059:            WRITE(MYUNIT2,'(3F28.20)') QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)
1061:         ENDDO1060:         ENDDO
1062: 1061:         
1063:         CLOSE(MYUNIT2)1062:         CLOSE(MYUNIT2)
1064: 1063:         
1065:      ENDDO1064:      ENDDO
1066:      !1065:      !
1067:      !  Write out lowest NSAVE structures to xmakemol ellipsoid format for1066:      !  Write out lowest NSAVE structures to xmakemol ellipsoid format for
1068:      !  clusters of ellipsoids of revolution1067:      !  clusters of ellipsoids of revolution
1069:      !1068:      !
1070: 1069:      
1071:      IF (MPIT.OR.GBHT) THEN1070:      IF (MPIT.OR.GBHT) THEN
1072:         WRITE (ISTR, '(I10)') MYNODE+11071:         WRITE (ISTR, '(I10)') MYNODE+1
1073:         MYUNIT2=GETUNIT()1072:         MYUNIT2=GETUNIT()
1074:         MYFILENAME2="ellipsoid."//trim(adjustl(istr))//".xyz"1073:         MYFILENAME2="ellipsoid."//trim(adjustl(istr))//".xyz"
1075:         OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")1074:         OPEN(MYUNIT2,FILE=trim(adjustl(MYFILENAME2)), STATUS="unknown", form="formatted")
1076:      ELSE1075:      ELSE
1077:         MYUNIT2=GETUNIT()1076:         MYUNIT2=GETUNIT()
1078:         OPEN(MYUNIT2,FILE='ellipsoid.xyz',STATUS='UNKNOWN')1077:         OPEN(MYUNIT2,FILE='ellipsoid.xyz',STATUS='UNKNOWN')
1079:      ENDIF1078:      ENDIF
1080: 1079:      
1081:      IF(PYT) THEN1080:      IF(PYT) THEN
1082:         CALL PY_OUTPUT(MYUNIT2)1081:         CALL PY_OUTPUT(MYUNIT2)
1083:      ELSE1082:      ELSE
1084: 1083:         
1085:         do J1=1,NSAVE1084:         do J1=1,NSAVE
1086:            WRITE(MYUNIT2,*) NATOMS/21085:            WRITE(MYUNIT2,*) NATOMS/2
1087:            WRITE(MYUNIT2,10) J1, QMIN(J1), FF(J1)1086:            WRITE(MYUNIT2,10) J1, QMIN(J1), FF(J1)
1088: 1087:            
1089:            DO J2=1,NATOMS/21088:            DO J2=1,NATOMS/2
1090: 1089:               
1091:               IF (GAYBERNET) THEN1090:               IF (GAYBERNET) THEN
1092:                  CALL EllipsoidsAAtoPolar(QMINP(J1,3*NATOMS/2+3*(J2-1)+1),QMINP(J1,3*NATOMS/2+3*(J2-1)+2),&1091:                  CALL EllipsoidsAAtoPolar(QMINP(J1,3*NATOMS/2+3*(J2-1)+1),QMINP(J1,3*NATOMS/2+3*(J2-1)+2),&
1093:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+3),EulerPhi,EulerPsi,EulerTheta,EulerPhiDeg,EulerPsiDeg,EulerThetaDeg)1092:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+3),EulerPhi,EulerPsi,EulerTheta,EulerPhiDeg,EulerPsiDeg,EulerThetaDeg)
1094:                  !  EulerPhiDeg = 90-EulerPhiDeg    ! EulerPhiDeg returned from EllipsoidsAAtoPolar is in fact the alpha angle1093:                  !  EulerPhiDeg = 90-EulerPhiDeg    ! EulerPhiDeg returned from EllipsoidsAAtoPolar is in fact the alpha angle
1095:                  !                                  ! defined by me (angle of vector with the xy plane)1094:                  !                                  ! defined by me (angle of vector with the xy plane)
1096:                  WRITE(MYUNIT2,'(a5,2x,3f20.10,2x,a8,6f15.8,2x,a11,3f15.8)') '0',QMINP(J1,3*(J2-1)+1),&1095:                  WRITE(MYUNIT2,'(a5,2x,3f20.10,2x,a8,6f15.8,2x,a11,3f15.8)') '0',QMINP(J1,3*(J2-1)+1),&
1097:                       QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&1096:                       QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&
1098:                       'ellipse ',1.0D0,1.0D0,GBANISOTROPYR,&1097:                       'ellipse ',1.0D0,1.0D0,GBANISOTROPYR,&
1099:                       EulerPsiDeg,EulerPhiDeg,0.0D0,&! this is in degrees1098:                       EulerPsiDeg,EulerPhiDeg,0.0D0,&! this is in degrees
1100:                       'atom_vector',QMINP(J1,3*NATOMS/2+3*(J2-1)+1),&1099:                       'atom_vector',QMINP(J1,3*NATOMS/2+3*(J2-1)+1),&
1101:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+2),QMINP(J1,3*NATOMS/2+3*(J2-1)+3)1100:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+2),QMINP(J1,3*NATOMS/2+3*(J2-1)+3)
1102:               ELSE IF (PYGPERIODICT.OR.PYBINARYT) THEN1101:               ELSE IF (PYGPERIODICT.OR.PYBINARYT) THEN
1103:                  CALL AAtoEuler(QMINP(J1,3*NATOMS/2+3*(J2-1)+1),QMINP(J1,3*NATOMS/2+3*(J2-1)+2),&1102:                  CALL AAtoEuler(QMINP(J1,3*NATOMS/2+3*(J2-1)+1),QMINP(J1,3*NATOMS/2+3*(J2-1)+2),&
1104:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+3),EulerPhiDeg,EulerPsiDeg,EulerThetaDeg)1103:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+3),EulerPhiDeg,EulerPsiDeg,EulerThetaDeg)
1105: 1104:                  
1106:                  WRITE(MYUNIT2,'(a5,2x,3f20.10,2x,a8,6f15.8,2x,a11,3f15.8)') 'O',QMINP(J1,3*(J2-1)+1),&1105:                  WRITE(MYUNIT2,'(a5,2x,3f20.10,2x,a8,6f15.8,2x,a11,3f15.8)') 'O',QMINP(J1,3*(J2-1)+1),&
1107:                       QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&1106:                       QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&
1108:                       'ellipse ',PYA1BIN(J2,1)*2.0D0,PYA1BIN(J2,2)*2.0D0,PYA1BIN(J2,3)*2.0D0,EulerPhiDeg,EulerPsiDeg,EulerThetaDeg,&1107:                       'ellipse ',PYA1BIN(J2,1)*2.0D0,PYA1BIN(J2,2)*2.0D0,PYA1BIN(J2,3)*2.0D0,EulerPhiDeg,EulerPsiDeg,EulerThetaDeg,&
1109:                       'atom_vector',QMINP(J1,3*NATOMS/2+3*(J2-1)+1),&1108:                       'atom_vector',QMINP(J1,3*NATOMS/2+3*(J2-1)+1),&
1110:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+2),QMINP(J1,3*NATOMS/2+3*(J2-1)+3)1109:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+2),QMINP(J1,3*NATOMS/2+3*(J2-1)+3)
1111:               ELSE IF (GBT.OR.GBDT) THEN1110:               ELSE IF (GBT.OR.GBDT) THEN
1112:                  CALL AAtoEuler(QMINP(J1,3*NATOMS/2+3*(J2-1)+1),QMINP(J1,3*NATOMS/2+3*(J2-1)+2),&1111:                  CALL AAtoEuler(QMINP(J1,3*NATOMS/2+3*(J2-1)+1),QMINP(J1,3*NATOMS/2+3*(J2-1)+2),&
1113:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+3),EulerPhiDeg,EulerPsiDeg,EulerThetaDeg)1112:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+3),EulerPhiDeg,EulerPsiDeg,EulerThetaDeg)
1114: 1113:                  
1115:                  WRITE(MYUNIT2,'(a5,2x,3f20.10,2x,a8,6f15.8,2x,a11,3f15.8)') 'O',QMINP(J1,3*(J2-1)+1),&1114:                  WRITE(MYUNIT2,'(a5,2x,3f20.10,2x,a8,6f15.8,2x,a11,3f15.8)') 'O',QMINP(J1,3*(J2-1)+1),&
1116:                       QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&1115:                       QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&
1117:                       'ellipse ',GBKAPPA,1.0D0,1.0D0,EulerPhiDeg,EulerPsiDeg,EulerThetaDeg,&1116:                       'ellipse ',GBKAPPA,1.0D0,1.0D0,EulerPhiDeg,EulerPsiDeg,EulerThetaDeg,&
1118:                       'atom_vector',QMINP(J1,3*NATOMS/2+3*(J2-1)+1),&1117:                       'atom_vector',QMINP(J1,3*NATOMS/2+3*(J2-1)+1),&
1119:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+2),QMINP(J1,3*NATOMS/2+3*(J2-1)+3)1118:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+2),QMINP(J1,3*NATOMS/2+3*(J2-1)+3)
1120:               ELSE IF (LJCAPSIDT) THEN1119:               ELSE IF (LJCAPSIDT) THEN
1121:                  CALL AAtoEuler(QMINP(J1,3*NATOMS/2+3*(J2-1)+1),QMINP(J1,3*NATOMS/2+3*(J2-1)+2),&1120:                  CALL AAtoEuler(QMINP(J1,3*NATOMS/2+3*(J2-1)+1),QMINP(J1,3*NATOMS/2+3*(J2-1)+2),&
1122:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+3),EulerPhiDeg,EulerPsiDeg,EulerThetaDeg)1121:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+3),EulerPhiDeg,EulerPsiDeg,EulerThetaDeg)
1123: 1122:                  
1124:                  WRITE(MYUNIT2,'(a5,2x,3f20.10,2x,a8,6f15.8,2x,a11,3f15.8)') 'O',QMINP(J1,3*(J2-1)+1),&1123:                  WRITE(MYUNIT2,'(a5,2x,3f20.10,2x,a8,6f15.8,2x,a11,3f15.8)') 'O',QMINP(J1,3*(J2-1)+1),&
1125:                       QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&1124:                       QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&
1126:                       'ellipse ',1.0D0,1.0D0,1.0D0,EulerPhiDeg,EulerPsiDeg,EulerThetaDeg,&1125:                       'ellipse ',1.0D0,1.0D0,1.0D0,EulerPhiDeg,EulerPsiDeg,EulerThetaDeg,&
1127:                       'atom_vector',QMINP(J1,3*NATOMS/2+3*(J2-1)+1),&1126:                       'atom_vector',QMINP(J1,3*NATOMS/2+3*(J2-1)+1),&
1128:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+2),QMINP(J1,3*NATOMS/2+3*(J2-1)+3)1127:                       QMINP(J1,3*NATOMS/2+3*(J2-1)+2),QMINP(J1,3*NATOMS/2+3*(J2-1)+3)
1129:               ENDIF1128:               ENDIF
1130:            ENDDO1129:            ENDDO
1131:         ENDDO1130:         ENDDO
1132: 1131:         
1133:      ENDIF ! PYT1132:      ENDIF ! PYT
1134: 1133:      
1135:      CLOSE(MYUNIT2)1134:      CLOSE(MYUNIT2)
1136: 1135:      
1137:   ELSE IF (TIP) THEN1136:   ELSE IF (TIP) THEN
1138:      LUNIT=GETUNIT()1137:      LUNIT=GETUNIT()
1139:      OPEN(UNIT=LUNIT,FILE='tip.xyz',STATUS='UNKNOWN')1138:      OPEN(UNIT=LUNIT,FILE='tip.xyz',STATUS='UNKNOWN')
1140:      DO J1=1,NSAVE1139:      DO J1=1,NSAVE
1141:         WRITE(LUNIT,'(I6)') (NATOMS/2)*31140:         WRITE(LUNIT,'(I6)') (NATOMS/2)*3
1142:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1141:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1143:         DO J2=1,NATOMS/21142:         DO J2=1,NATOMS/2
1144:            CALL TIPIO(QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&1143:            CALL TIPIO(QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&
1145:                 QMINP(J1,3*(NATOMS/2+J2-1)+1),QMINP(J1,3*(NATOMS/2+J2-1)+2),QMINP(J1,3*(NATOMS/2+J2-1)+3),&1144:                 QMINP(J1,3*(NATOMS/2+J2-1)+1),QMINP(J1,3*(NATOMS/2+J2-1)+2),QMINP(J1,3*(NATOMS/2+J2-1)+3),&
1146:                 RBCOORDS)1145:                 RBCOORDS)
1198:                     ELSE1197:                     ELSE
1199:                        WRITE (LUNIT,'(A4,3F18.10,A12,3F18.10)') 'LB ',QMINP(J1,3*(J4-1)+1),QMINP(J1,3*(J4-1)+2),QMINP(J1,3*(J4-1)+3)1198:                        WRITE (LUNIT,'(A4,3F18.10,A12,3F18.10)') 'LB ',QMINP(J1,3*(J4-1)+1),QMINP(J1,3*(J4-1)+2),QMINP(J1,3*(J4-1)+3)
1200:                     ENDIF1199:                     ENDIF
1201:                  ENDDO1200:                  ENDDO
1202:               ENDIF1201:               ENDIF
1203:               NDUMMY=NDUMMY+NSITEPERBODY(J2)1202:               NDUMMY=NDUMMY+NSITEPERBODY(J2)
1204:            ENDDO1203:            ENDDO
1205:         ENDDO1204:         ENDDO
1206:         CLOSE(LUNIT)1205:         CLOSE(LUNIT)
1207:      ENDIF1206:      ENDIF
1208: 1207:      
1209:   ELSE IF (DBPT .OR. DMBLPYT) THEN1208:   ELSE IF (DBPT .OR. DMBLPYT) THEN
1210: 1209:      
1211:      LUNIT=GETUNIT()1210:      LUNIT=GETUNIT()
1212:      IF (DBPT) OPEN(UNIT=LUNIT, FILE='dbp.xyz', STATUS='UNKNOWN')1211:      IF (DBPT) OPEN(UNIT=LUNIT, FILE='dbp.xyz', STATUS='UNKNOWN')
1213:      IF (DMBLPYT) OPEN(UNIT=LUNIT, FILE='dmblpy.xyz', STATUS='UNKNOWN')1212:      IF (DMBLPYT) OPEN(UNIT=LUNIT, FILE='dmblpy.xyz', STATUS='UNKNOWN')
1214: 1213:      
1215:      GTEST = .FALSE.1214:      GTEST = .FALSE.
1216:      DU    = (/0.D0, 1.D0, 0.D0/)1215:      DU    = (/0.D0, 1.D0, 0.D0/)
1217:      DO J1 = 1, NSAVE1216:      DO J1 = 1, NSAVE
1218: 1217:                 
1219:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES1218:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES 
1220:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1219:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1221: 1220:         
1222:         DO J2 = 1, NATOMS/21221:         DO J2 = 1, NATOMS/2
1223: 1222:            
1224:            J3   = 3*J21223:            J3   = 3*J2
1225:            J5   = 3*NATOMS/2 + J31224:            J5   = 3*NATOMS/2 + J3
1226:            P(:) = QMINP(J1,J5-2:J5)1225:            P(:) = QMINP(J1,J5-2:J5)
1227: 1226:            
1228:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)1227:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)
1229: 1228:            
1230:            DO J4 = 1, NRBSITES1229:            DO J4 = 1, NRBSITES
1231: 1230:               
1232:               IF (J4 == 1) THEN1231:               IF (J4 == 1) THEN
1233:                  RBCOORDS(1:3) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))1232:                  RBCOORDS(1:3) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))
1234:                  WRITE(LUNIT,'(A4,3F20.10)') 'O', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)1233:                  WRITE(LUNIT,'(A4,3F20.10)') 'O', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)
1235:               ELSEIF (J4 == 2) THEN1234:               ELSEIF (J4 == 2) THEN
1236:                  RBCOORDS(1:3) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))1235:                  RBCOORDS(1:3) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))
1237:                  WRITE(LUNIT,'(A4,3F20.10)') 'C', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)1236:                  WRITE(LUNIT,'(A4,3F20.10)') 'C', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)
1238:               ELSE1237:               ELSE
1239:                  RBCOORDS(1:3) = MATMUL(RMI(:,:),DU)1238:                  RBCOORDS(1:3) = MATMUL(RMI(:,:),DU)
1240:                  WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&1239:                  WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&
1241:                       'H', QMINP(J1,J3-2), QMINP(J1,J3-1), QMINP(J1,J3),&1240:                       'H', QMINP(J1,J3-2), QMINP(J1,J3-1), QMINP(J1,J3),&
1242:                       'atom_vector', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)1241:                       'atom_vector', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)
1243:               ENDIF1242:               ENDIF
1244: 1243:               
1245:            ENDDO1244:            ENDDO
1246: 1245:            
1247:         ENDDO1246:         ENDDO
1248: 1247:         
1249:      ENDDO1248:      ENDDO
1250:      CLOSE(LUNIT)1249:      CLOSE(LUNIT)
1251: 1250:      
1252:      RETURN1251:      RETURN
1253: 1252:      
1254:   ELSE IF (DBPTDT) THEN1253:   ELSE IF (DBPTDT) THEN
1255: 1254:      
1256:      CALL VIEWDMBLTD()1255:      CALL VIEWDMBLTD()
1257:      RETURN1256:      RETURN
1258: 1257:      
1259:   ELSE IF (DMBLMT) THEN1258:   ELSE IF (DMBLMT) THEN
1260: 1259:      
1261:      CALL VIEWDMBL()1260:      CALL VIEWDMBL()
1262:      RETURN1261:      RETURN
1263: 1262:      
1264:   ELSE IF (GBT .OR. GBDT) THEN1263:   ELSE IF (GBT .OR. GBDT) THEN
1265: 1264:      
1266:      LUNIT=GETUNIT()1265:      LUNIT=GETUNIT()
1267:      OPEN(UNIT=LUNIT, FILE='gbe.xyz', STATUS='UNKNOWN')1266:      OPEN(UNIT=LUNIT, FILE='gbe.xyz', STATUS='UNKNOWN')
1268:      GTEST = .FALSE.1267:      GTEST = .FALSE.
1269:      DO J1 = 1, NSAVE1268:      DO J1 = 1, NSAVE
1270: 1269:         
1271:         WRITE(LUNIT,'(I6)') NATOMS/21270:         WRITE(LUNIT,'(I6)') NATOMS/2
1272:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1271:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1273: 1272:         
1274:         DO J2 = 1, NATOMS/21273:         DO J2 = 1, NATOMS/2
1275: 1274:            
1276:            J3   = 3*J21275:            J3   = 3*J2
1277:            J5   = 3*NATOMS/2 + J31276:            J5   = 3*NATOMS/2 + J3
1278:            RBCOORDS(1:3) = QMINP(J1,J3-2:J3)1277:            RBCOORDS(1:3) = QMINP(J1,J3-2:J3)
1279:            P(:)          = QMINP(J1,J5-2:J5)1278:            P(:)          = QMINP(J1,J5-2:J5)
1280: 1279:            
1281:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)1280:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)
1282: 1281:            
1283:            PHI   = DATAN2(RMI(2,3),RMI(1,3))1282:            PHI   = DATAN2(RMI(2,3),RMI(1,3))
1284:            IF (PHI <= 0.D0) PHI = PHI + 2.D0*PI1283:            IF (PHI <= 0.D0) PHI = PHI + 2.D0*PI
1285: 1284:            
1286:            THT   = DACOS(RMI(3,3))1285:            THT   = DACOS(RMI(3,3))
1287: 1286:            
1288:            PHI   = PHI*180.D0/PI1287:            PHI   = PHI*180.D0/PI
1289:            THT   = THT*180.D0/PI1288:            THT   = THT*180.D0/PI
1290: 1289:            
1291:            WRITE(LUNIT,'(a5,2x,3f20.10,2x,a8,6f20.10)')&1290:            WRITE(LUNIT,'(a5,2x,3f20.10,2x,a8,6f20.10)')&
1292:                 'O', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3),&1291:                 'O', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3),&
1293:                 'ellipse', 2.D0*ESA(1), 2.D0*ESA(2), 2.D0*ESA(3), PHI, THT, 0.D01292:                 'ellipse', 2.D0*ESA(1), 2.D0*ESA(2), 2.D0*ESA(3), PHI, THT, 0.D0
1294: 1293:            
1295:         ENDDO1294:         ENDDO
1296: 1295:         
1297:      ENDDO1296:      ENDDO
1298: 1297:      
1299:      RETURN1298:      RETURN
1300: 1299:      
1301:   ELSE IF (LINRODT) THEN1300:   ELSE IF (LINRODT) THEN
1302: 1301:      
1303:      LUNIT=GETUNIT()1302:      LUNIT=GETUNIT()
1304:      OPEN(UNIT=LUNIT, FILE='linrod.xyz', STATUS='UNKNOWN')1303:      OPEN(UNIT=LUNIT, FILE='linrod.xyz', STATUS='UNKNOWN')
1305:      GTEST = .FALSE.1304:      GTEST = .FALSE.
1306: 1305:      
1307:      DO J1 = 1, NSAVE1306:      DO J1 = 1, NSAVE
1308: 1307:         
1309:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES1308:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES
1310:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1309:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1311: 1310:         
1312:         DO J2 = 1, NATOMS/21311:         DO J2 = 1, NATOMS/2
1313: 1312:            
1314:            J3    = 3*J21313:            J3    = 3*J2
1315:            J5    = 3*NATOMS/2 + J31314:            J5    = 3*NATOMS/2 + J3
1316:            P(:)  = QMINP(J1,J5-2:J5)1315:            P(:)  = QMINP(J1,J5-2:J5)
1317: 1316:            
1318:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)1317:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)
1319: 1318:            
1320: 1319:            
1321:            DO J4 = 1, NRBSITES1320:            DO J4 = 1, NRBSITES
1322: 1321:               
1323:               RBCOORDS(3*J4-2:3*J4) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))1322:               RBCOORDS(3*J4-2:3*J4) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))
1324: 1323:               
1325:            ENDDO1324:            ENDDO
1326: 1325:            
1327:            DO J4 = 1, NRBSITES1326:            DO J4 = 1, NRBSITES
1328: 1327:               
1329:               J3 = J4 + 11328:               J3 = J4 + 1
1330:               IF (J4 == NRBSITES) J3 = 11329:               IF (J4 == NRBSITES) J3 = 1
1331:               P(:) = RBCOORDS(3*J3-2:3*J3) - RBCOORDS(3*J4-2:3*J4)1330:               P(:) = RBCOORDS(3*J3-2:3*J3) - RBCOORDS(3*J4-2:3*J4)
1332:               WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&1331:               WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&
1333:                    'O', RBCOORDS(3*J4-2), RBCOORDS(3*J4-1), RBCOORDS(3*J4), 'atom_vector', P(1), P(2), P(3)1332:                    'O', RBCOORDS(3*J4-2), RBCOORDS(3*J4-1), RBCOORDS(3*J4), 'atom_vector', P(1), P(2), P(3)
1334: 1333:               
1335:            ENDDO1334:            ENDDO
1336: 1335:            
1337:         ENDDO1336:         ENDDO
1338: 1337:         
1339:      ENDDO1338:      ENDDO
1340: 1339:      
1341:      CLOSE(UNIT=LUNIT)1340:      CLOSE(UNIT=LUNIT)
1342: 1341:      
1343:      RETURN1342:      RETURN
1344: 1343:      
1345:   ELSE IF (LWOTPT) THEN1344:   ELSE IF (LWOTPT) THEN
1346: 1345:      
1347:      LUNIT=GETUNIT()1346:      LUNIT=GETUNIT()
1348:      OPEN(UNIT=LUNIT, FILE='lwotp.xyz', STATUS='UNKNOWN')1347:      OPEN(UNIT=LUNIT, FILE='lwotp.xyz', STATUS='UNKNOWN')
1349:      GTEST = .FALSE.1348:      GTEST = .FALSE.
1350: 1349:      
1351:      DO J1 = 1, NSAVE1350:      DO J1 = 1, NSAVE
1352: 1351:         
1353:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES1352:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES
1354:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1353:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1355: 1354:         
1356:         DO J2 = 1, NATOMS/21355:         DO J2 = 1, NATOMS/2
1357: 1356:            
1358:            J3    = 3*J21357:            J3    = 3*J2
1359:            J5    = 3*NATOMS/2 + J31358:            J5    = 3*NATOMS/2 + J3
1360:            P(:)  = QMINP(J1,J5-2:J5)1359:            P(:)  = QMINP(J1,J5-2:J5)
1361: 1360:            
1362:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)1361:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)
1363: 1362:            
1364: 1363:            
1365:            DO J4 = 1, NRBSITES1364:            DO J4 = 1, NRBSITES
1366: 1365:               
1367:               RBCOORDS(3*J4-2:3*J4) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))1366:               RBCOORDS(3*J4-2:3*J4) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))
1368: 1367:               
1369:            ENDDO1368:            ENDDO
1370: 1369:            
1371:            DO J4 = 1, NRBSITES1370:            DO J4 = 1, NRBSITES
1372: 1371:               
1373:               J3 = J4 + 11372:               J3 = J4 + 1
1374:               IF (J4 == NRBSITES) J3 = 11373:               IF (J4 == NRBSITES) J3 = 1
1375:               P(:) = RBCOORDS(3*J3-2:3*J3) - RBCOORDS(3*J4-2:3*J4)1374:               P(:) = RBCOORDS(3*J3-2:3*J3) - RBCOORDS(3*J4-2:3*J4)
1376:               WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&1375:               WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&
1377:                    'O', RBCOORDS(3*J4-2), RBCOORDS(3*J4-1), RBCOORDS(3*J4), 'atom_vector', P(1), P(2), P(3)1376:                    'O', RBCOORDS(3*J4-2), RBCOORDS(3*J4-1), RBCOORDS(3*J4), 'atom_vector', P(1), P(2), P(3)
1378: 1377:               
1379:            ENDDO1378:            ENDDO
1380: 1379:            
1381:         ENDDO1380:         ENDDO
1382: 1381:         
1383:      ENDDO1382:      ENDDO
1384: 1383:      
1385:      CLOSE(UNIT=LUNIT)1384:      CLOSE(UNIT=LUNIT)
1386: 1385:      
1387:      RETURN1386:      RETURN
1388: 1387:      
1389:   ELSE IF (MSTBINT) THEN1388:   ELSE IF (MSTBINT) THEN
1390: 1389:      
1391:      LUNIT=GETUNIT()1390:      LUNIT=GETUNIT()
1392:      OPEN(UNIT=LUNIT, FILE='mstbin.xyz', STATUS='UNKNOWN')1391:      OPEN(UNIT=LUNIT, FILE='mstbin.xyz', STATUS='UNKNOWN')
1393:      GTEST = .FALSE.1392:      GTEST = .FALSE.
1394: 1393:      
1395:      DO J1 = 1, NSAVE1394:      DO J1 = 1, NSAVE
1396: 1395:         
1397:         WRITE(LUNIT,'(I6)') NPS*NRBSITES1 + (NATOMS/2 - NPS)*(NRBSITES - NRBSITES1)1396:         WRITE(LUNIT,'(I6)') NPS*NRBSITES1 + (NATOMS/2 - NPS)*(NRBSITES - NRBSITES1)
1398:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1397:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1399: 1398:         
1400:         DO J2 = 1, NATOMS/21399:         DO J2 = 1, NATOMS/2
1401: 1400:            
1402:            J3    = 3*J21401:            J3    = 3*J2
1403:            J5    = 3*NATOMS/2 + J31402:            J5    = 3*NATOMS/2 + J3
1404:            P(:)  = QMINP(J1,J5-2:J5)1403:            P(:)  = QMINP(J1,J5-2:J5)
1405: 1404:            
1406:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)1405:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)
1407: 1406:            
1408:            IF (J2 <= NPS) THEN1407:            IF (J2 <= NPS) THEN
1409:               NRBS1 = 11408:               NRBS1 = 1
1410:               NRBS2 = NRBSITES11409:               NRBS2 = NRBSITES1
1411:            ELSE1410:            ELSE
1412:               NRBS1 = NRBSITES1 + 11411:               NRBS1 = NRBSITES1 + 1
1413:               NRBS2 = NRBSITES1412:               NRBS2 = NRBSITES
1414:            ENDIF1413:            ENDIF
1415: 1414:            
1416:            DO J4 = NRBS1, NRBS21415:            DO J4 = NRBS1, NRBS2
1417: 1416:               
1418:               RBCOORDS(3*J4-2:3*J4) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))1417:               RBCOORDS(3*J4-2:3*J4) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))
1419: 1418:               
1420:            ENDDO1419:            ENDDO
1421: 1420:            
1422:            DO J4 = NRBS1, NRBS21421:            DO J4 = NRBS1, NRBS2
1423: 1422:               
1424:               J3 = J4 + 11423:               J3 = J4 + 1
1425: 1424:               
1426:               IF (J2 <= NPS) THEN1425:               IF (J2 <= NPS) THEN
1427:                  IF (J4 == NRBSITES1) J3 = 11426:                  IF (J4 == NRBSITES1) J3 = 1
1428:               ELSE1427:               ELSE
1429:                  IF (J4 == NRBSITES) J3 = NRBSITES1 + 11428:                  IF (J4 == NRBSITES) J3 = NRBSITES1 + 1
1430:               ENDIF1429:               ENDIF
1431: 1430:               
1432:               P(:) = RBCOORDS(3*J3-2:3*J3) - RBCOORDS(3*J4-2:3*J4)1431:               P(:) = RBCOORDS(3*J3-2:3*J3) - RBCOORDS(3*J4-2:3*J4)
1433:               IF (J2 <= NPS) THEN1432:               IF (J2 <= NPS) THEN
1434:                  WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&1433:                  WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&
1435:                       'O', RBCOORDS(3*J4-2), RBCOORDS(3*J4-1), RBCOORDS(3*J4), 'atom_vector', P(1), P(2), P(3)1434:                       'O', RBCOORDS(3*J4-2), RBCOORDS(3*J4-1), RBCOORDS(3*J4), 'atom_vector', P(1), P(2), P(3)
1436:               ELSE1435:               ELSE
1437:                  WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&1436:                  WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&
1438:                       'N', RBCOORDS(3*J4-2), RBCOORDS(3*J4-1), RBCOORDS(3*J4), 'atom_vector', P(1), P(2), P(3)1437:                       'N', RBCOORDS(3*J4-2), RBCOORDS(3*J4-1), RBCOORDS(3*J4), 'atom_vector', P(1), P(2), P(3)
1439:               ENDIF1438:               ENDIF
1440: 1439:               
1441:            ENDDO1440:            ENDDO
1442: 1441:            
1443:         ENDDO1442:         ENDDO
1444: 1443:         
1445:      ENDDO1444:      ENDDO
1446: 1445:      
1447:      CLOSE(UNIT=LUNIT)1446:      CLOSE(UNIT=LUNIT)
1448: 1447:      
1449:      RETURN1448:      RETURN
1450: 1449:      
1451:   ELSE IF (MSSTOCKT) THEN1450:   ELSE IF (MSSTOCKT) THEN
1452: 1451:      
1453:      LUNIT=GETUNIT()1452:      LUNIT=GETUNIT()
1454:      OPEN(UNIT=LUNIT, FILE='msstock.xyz', STATUS='UNKNOWN')1453:      OPEN(UNIT=LUNIT, FILE='msstock.xyz', STATUS='UNKNOWN')
1455:      GTEST = .FALSE.1454:      GTEST = .FALSE.
1456: 1455:      
1457:      DO J1 = 1, NSAVE1456:      DO J1 = 1, NSAVE
1458: 1457:         
1459:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES1458:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES
1460:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1459:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1461: 1460:         
1462:         DO J2 = 1, NATOMS/21461:         DO J2 = 1, NATOMS/2
1463: 1462:            
1464:            J3   = 3*J21463:            J3   = 3*J2
1465:            J5   = 3*NATOMS/2 + J31464:            J5   = 3*NATOMS/2 + J3
1466:            P(:) = QMINP(J1,J5-2:J5)1465:            P(:) = QMINP(J1,J5-2:J5)
1467: 1466:            
1468:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)1467:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)
1469: 1468:            
1470:            DO J4 = 1, NRBSITES1469:            DO J4 = 1, NRBSITES
1471: 1470:               
1472:               RBCOORDS(1:3) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))1471:               RBCOORDS(1:3) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))
1473:               P(:)          = MATMUL(RMI(:,:),RBUV(J4,:))1472:               P(:)          = MATMUL(RMI(:,:),RBUV(J4,:))
1474:               WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&1473:               WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&
1475:                    'O', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3), 'atom_vector', P(1), P(2), P(3)1474:                    'O', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3), 'atom_vector', P(1), P(2), P(3)
1476: 1475:               
1477:            ENDDO1476:            ENDDO
1478: 1477:            
1479:         ENDDO1478:         ENDDO
1480: 1479:         
1481:      ENDDO1480:      ENDDO
1482: 1481:      
1483:      CLOSE(UNIT=LUNIT)1482:      CLOSE(UNIT=LUNIT)
1484: 1483:      
1485:      LUNIT=GETUNIT()1484:      LUNIT=GETUNIT()
1486:      OPEN(UNIT=LUNIT, FILE='msstktr.xyz', STATUS='UNKNOWN')1485:      OPEN(UNIT=LUNIT, FILE='msstktr.xyz', STATUS='UNKNOWN')
1487:      GTEST = .FALSE.1486:      GTEST = .FALSE.
1488: 1487:      
1489:      DO J1 = 1, NSAVE1488:      DO J1 = 1, NSAVE
1490: 1489:         
1491:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES !(NRBSITES - 1)1490:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES !(NRBSITES - 1)
1492:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1491:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1493: 1492:         
1494:         DO J2 = 1, NATOMS/21493:         DO J2 = 1, NATOMS/2
1495: 1494:            
1496:            J3    = 3*J21495:            J3    = 3*J2
1497:            J5    = 3*NATOMS/2 + J31496:            J5    = 3*NATOMS/2 + J3
1498:            P(:)  = QMINP(J1,J5-2:J5)1497:            P(:)  = QMINP(J1,J5-2:J5)
1499: 1498:            
1500:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)1499:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)
1501: 1500:            
1502:            DO J4 = 1, NRBSITES1501:            DO J4 = 1, NRBSITES
1503: 1502:               
1504:               RBCOORDS(3*J4-2:3*J4) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))1503:               RBCOORDS(3*J4-2:3*J4) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))
1505: 1504:               
1506:            ENDDO1505:            ENDDO
1507: 1506:            
1508:            DO J4 = 1, NRBSITES !- 11507:            DO J4 = 1, NRBSITES !- 1
1509: 1508:               
1510:               J3 = J4 + 11509:               J3 = J4 + 1
1511:               IF (J4 == NRBSITES) J3 = 11510:               IF (J4 == NRBSITES) J3 = 1
1512:               !                  IF (J4 == NRBSITES - 1) J3 = 11511:               !                  IF (J4 == NRBSITES - 1) J3 = 1
1513:               P(:) = RBCOORDS(3*J3-2:3*J3) - RBCOORDS(3*J4-2:3*J4)1512:               P(:) = RBCOORDS(3*J3-2:3*J3) - RBCOORDS(3*J4-2:3*J4)
1514:               WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&1513:               WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)')&
1515:                    'O', RBCOORDS(3*J4-2), RBCOORDS(3*J4-1), RBCOORDS(3*J4), 'atom_vector', P(1), P(2), P(3)1514:                    'O', RBCOORDS(3*J4-2), RBCOORDS(3*J4-1), RBCOORDS(3*J4), 'atom_vector', P(1), P(2), P(3)
1516: 1515:               
1517:            ENDDO1516:            ENDDO
1518: 1517:            
1519:         ENDDO1518:         ENDDO
1520: 1519:         
1521:      ENDDO1520:      ENDDO
1522: 1521:      
1523:      CLOSE(UNIT=LUNIT)1522:      CLOSE(UNIT=LUNIT)
1524: 1523:      
1525:      RETURN1524:      RETURN
1526: 1525:      
1527:   ELSE IF (MULTPAHAT) THEN1526:   ELSE IF (MULTPAHAT) THEN
1528: 1527:      
1529:      CALL VIEWMULTPAHA()1528:      CALL VIEWMULTPAHA()
1530:      RETURN1529:      RETURN
1531: 1530:      
1532:   ELSE IF (NPAHT .OR. PAHAT .OR. PAHW99T) THEN1531:   ELSE IF (NPAHT .OR. PAHAT .OR. PAHW99T) THEN
1533: 1532:      
1534:      LUNIT=GETUNIT()1533:      LUNIT=GETUNIT()
1535:      OPEN(UNIT=LUNIT, FILE='rigid.xyz', STATUS='UNKNOWN')1534:      OPEN(UNIT=LUNIT, FILE='rigid.xyz', STATUS='UNKNOWN')
1536:      GTEST = .FALSE.1535:      GTEST = .FALSE.
1537: 1536:      
1538:      IF (PAHW99T) THEN1537:      IF (PAHW99T) THEN
1539:         NCPHST = NCARBON + (NRBSITES-NCARBON)/21538:         NCPHST = NCARBON + (NRBSITES-NCARBON)/2
1540:         DO J1 = 1, (NCPHST-NCARBON)1539:         DO J1 = 1, (NCPHST-NCARBON) 
1541:            SITE(NCARBON+J1,:) = SITE(NCPHST+J1,:)1540:            SITE(NCARBON+J1,:) = SITE(NCPHST+J1,:)
1542:         ENDDO1541:         ENDDO
1543:      ELSE1542:      ELSE
1544:         NCPHST = NRBSITES1543:         NCPHST = NRBSITES
1545:      ENDIF1544:      ENDIF
1546: 1545:      
1547:      DO J1 = 1, NSAVE1546:      DO J1 = 1, NSAVE
1548: 1547:         
1549:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES1548:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES
1550:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1549:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1551: 1550:         
1552:         DO J2 = 1, NATOMS/21551:         DO J2 = 1, NATOMS/2
1553: 1552:            
1554:            J3   = 3*J21553:            J3   = 3*J2
1555:            J5   = 3*NATOMS/2 + J31554:            J5   = 3*NATOMS/2 + J3
1556:            P(:) = QMINP(J1,J5-2:J5)1555:            P(:) = QMINP(J1,J5-2:J5)
1557: 1556:            
1558:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)1557:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)
1559: 1558:            
1560:            DO J4 = 1, NCPHST1559:            DO J4 = 1, NCPHST
1561: 1560:               
1562:               RBCOORDS(1:3) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))1561:               RBCOORDS(1:3) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))
1563:               IF (J4 <= NCARBON) THEN1562:               IF (J4 <= NCARBON) THEN
1564:                  WRITE(LUNIT,'(A4,3F20.10)') 'C', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)1563:                  WRITE(LUNIT,'(A4,3F20.10)') 'C', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)
1565:               ELSE1564:               ELSE
1566:                  WRITE(LUNIT,'(A4,3F20.10)') 'H', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)1565:                  WRITE(LUNIT,'(A4,3F20.10)') 'H', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)
1567:               ENDIF1566:               ENDIF
1568: 1567:               
1569:            ENDDO1568:            ENDDO
1570: 1569:            
1571:         ENDDO1570:         ENDDO
1572: 1571:         
1573:      ENDDO1572:      ENDDO
1574: 1573:      
1575:      RETURN1574:      RETURN
1576:   ELSE IF (CAPBINT) THEN1575:   ELSE IF (CAPBINT) THEN
1577:      CALL VIEWCAPBIN()1576:      CALL VIEWCAPBIN()
1578:      RETURN1577:      RETURN
1579: 1578:      
1580:   ELSE IF (NCAPT) THEN1579:   ELSE IF (NCAPT) THEN
1581: 1580:      
1582:      CALL VIEWNEWCAPSID()1581:      CALL VIEWNEWCAPSID()
1583:      RETURN1582:      RETURN
1584: 1583:      
1585:   ELSE IF (NTIPT) THEN1584:   ELSE IF (NTIPT) THEN
1586: 1585:      
1587:      CALL VIEWNEWTIP()1586:      CALL VIEWNEWTIP()
1588:      RETURN1587:      RETURN
1589: 1588:      
1590:   ELSE IF (MWFILMT) THEN1589:   ELSE IF (MWFILMT) THEN
1591: 1590:      
1592:      CALL MWDRAW(NATOMS, PERIODIC, LAT)1591:      CALL MWDRAW(NATOMS, PERIODIC, LAT)
1593:      RETURN1592:      RETURN
1594: 1593:      
1595:   ELSE IF (PAPT) THEN1594:   ELSE IF (PAPT) THEN
1596: 1595:      
1597:      CALL VIEWPAP()1596:      CALL VIEWPAP()
1598:      RETURN1597:      RETURN
1599: 1598:      
1600:   ELSE IF (PAPBINT) THEN1599:   ELSE IF (PAPBINT) THEN
1601: 1600:      
1602:      CALL VIEWPAPBIN()1601:      CALL VIEWPAPBIN()
1603:      RETURN1602:      RETURN
1604: 1603:      
1605:   ELSE IF (PAPJANT) THEN1604:   ELSE IF (PAPJANT) THEN
1606: 1605:      
1607:      CALL VIEWPAPJANUS()1606:      CALL VIEWPAPJANUS()
1608:      RETURN1607:      RETURN
1609: 1608: 
1610:   ELSE IF (POLYT) THEN1609:   ELSE IF (POLYT) THEN
1611: 1610: 
1612:      CALL VIEW_POLYHEDRA()1611:      CALL VIEW_POLYHEDRA()
1613:      RETURN1612:      RETURN
1614: 1613:      
1615:   ELSE IF (LJ_GAUSST) THEN 
1616:  
1617:      CALL VIEW_LJ_GAUSS() 
1618:      RETURN 
1619:  
1620:   ELSE IF (PTSTSTT) THEN1614:   ELSE IF (PTSTSTT) THEN
1621: 1615:      
1622:      CALL VIEWPTSTST()1616:      CALL VIEWPTSTST()
1623:      RETURN1617:      RETURN
1624: 1618:      
1625:      !|gd351>1619:      !|gd351>
1626: 1620:      
1627:   ELSE IF (PATCHY) THEN1621:   ELSE IF (PATCHY) THEN
1628: 1622:      
1629:      CALL VIEWPATCHY()1623:      CALL VIEWPATCHY()
1630:      RETURN1624:      RETURN
1631: 1625:      
1632:      !<gd351|1626:      !<gd351|
1633:   ELSE IF (SANDBOXT) THEN1627:   ELSE IF (SANDBOXT) THEN
1634: 1628:      
1635:      CALL SANDBOX_OUTPUT()1629:      CALL SANDBOX_OUTPUT()
1636: 1630:      
1637:   ELSE IF (SILANET) THEN1631:   ELSE IF (SILANET) THEN
1638: 1632:      
1639:      CALL VIEWSILANE()1633:      CALL VIEWSILANE()
1640:      RETURN1634:      RETURN
1641: 1635:      
1642:   ELSE IF (STOCKAAT .OR. MORSEDPT) THEN1636:   ELSE IF (STOCKAAT .OR. MORSEDPT) THEN
1643: 1637:      
1644:      LUNIT=GETUNIT()1638:      LUNIT=GETUNIT()
1645:      IF (STOCKAAT) OPEN(UNIT=LUNIT, FILE='stockaa.xyz', STATUS='UNKNOWN')1639:      IF (STOCKAAT) OPEN(UNIT=LUNIT, FILE='stockaa.xyz', STATUS='UNKNOWN')
1646:      IF (MORSEDPT) OPEN(UNIT=LUNIT, FILE='morsedp.xyz', STATUS='UNKNOWN')1640:      IF (MORSEDPT) OPEN(UNIT=LUNIT, FILE='morsedp.xyz', STATUS='UNKNOWN')
1647:      GTEST = .FALSE.1641:      GTEST = .FALSE.
1648:      DU    = (/0.D0, 0.D0, 1.D0/)1642:      DU    = (/0.D0, 0.D0, 1.D0/)
1649: 1643:      
1650:      DO J1 = 1, NSAVE1644:      DO J1 = 1, NSAVE
1651: 1645:         
1652:         WRITE(LUNIT,'(I6)') NATOMS/21646:         WRITE(LUNIT,'(I6)') NATOMS/2
1653:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1647:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1654: 1648:         
1655:         DO J2 = 1, NATOMS/21649:         DO J2 = 1, NATOMS/2
1656: 1650:                
1657:            J3   = 3*J21651:            J3   = 3*J2
1658:            J5   = 3*NATOMS/2 + J31652:            J5   = 3*NATOMS/2 + J3
1659:            P(:) = QMINP(J1,J5-2:J5)1653:            P(:) = QMINP(J1,J5-2:J5)
1660: 1654:            
1661:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)1655:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)
1662: 1656:            
1663:            RBCOORDS(1:3) = MATMUL(RMI(:,:),DU(:))1657:            RBCOORDS(1:3) = MATMUL(RMI(:,:),DU(:))
1664:            WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)') 'O', QMINP(J1,J3-2), QMINP(J1,J3-1), QMINP(J1,J3),&1658:            WRITE(LUNIT,'(A4,3F20.10,2X,A12,2X,3F20.10)') 'O', QMINP(J1,J3-2), QMINP(J1,J3-1), QMINP(J1,J3),&
1665:                     'atom_vector', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)1659:                     'atom_vector', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)
1666: 1660:            
1667: 1661:            
1668:         ENDDO1662:         ENDDO
1669: 1663:         
1670:      ENDDO1664:      ENDDO
1671: 1665:      
1672:      CLOSE (UNIT=LUNIT)1666:      CLOSE (UNIT=LUNIT)
1673: 1667:          
1674:      RETURN1668:      RETURN
1675: 1669:      
1676:   ELSE IF (TDHDT) THEN1670:   ELSE IF (TDHDT) THEN
1677: 1671:      
1678:      CALL VIEWTDHD()1672:      CALL VIEWTDHD()
1679:      RETURN1673:      RETURN
1680: 1674:      
1681:   ELSE IF (DDMT) THEN1675:   ELSE IF (DDMT) THEN
1682: 1676:          
1683:      CALL VIEWDDM()1677:      CALL VIEWDDM()
1684:      RETURN1678:      RETURN
1685: 1679:      
1686:   ELSE IF (WATERDCT .OR. WATERKZT) THEN1680:   ELSE IF (WATERDCT .OR. WATERKZT) THEN
1687: 1681:      
1688:      LUNIT=GETUNIT()1682:      LUNIT=GETUNIT()
1689:      OPEN(UNIT=LUNIT, FILE='rigid.xyz', STATUS='UNKNOWN')1683:      OPEN(UNIT=LUNIT, FILE='rigid.xyz', STATUS='UNKNOWN')
1690:      GTEST = .FALSE.1684:      GTEST = .FALSE.
1691: 1685:      
1692:      DO J1 = 1, NSAVE1686:      DO J1 = 1, NSAVE
1693: 1687:         
1694:         WRITE(LUNIT,'(I6)') (NATOMS/2)*(NRBSITES - 1)1688:         WRITE(LUNIT,'(I6)') (NATOMS/2)*(NRBSITES - 1)
1695:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1689:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1696: 1690:         
1697:         DO J2 = 1, NATOMS/21691:         DO J2 = 1, NATOMS/2
1698: 1692:            
1699:            J3   = 3*J21693:            J3   = 3*J2
1700:            J5   = 3*NATOMS/2 + J31694:            J5   = 3*NATOMS/2 + J3
1701:            P(:) = QMINP(J1,J5-2:J5)1695:            P(:) = QMINP(J1,J5-2:J5)
1702: 1696:            
1703:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)1697:            CALL RMDRVT(P, RMI, DRMI, DRMI, DRMI, GTEST)
1704: 1698:            
1705:            DO J4 = 1, NRBSITES - 11699:            DO J4 = 1, NRBSITES - 1
1706: 1700:                   
1707:               RBCOORDS(1:3) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))1701:               RBCOORDS(1:3) = QMINP(J1,J3-2:J3) + MATMUL(RMI(:,:),SITE(J4,:))
1708:               IF (J4 == 1) THEN1702:               IF (J4 == 1) THEN
1709:                  WRITE(LUNIT,'(A4,3F20.10)') 'O', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)1703:                  WRITE(LUNIT,'(A4,3F20.10)') 'O', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)
1710:               ELSE1704:               ELSE
1711:                  WRITE(LUNIT,'(A4,3F20.10)') 'H', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)1705:                  WRITE(LUNIT,'(A4,3F20.10)') 'H', RBCOORDS(1), RBCOORDS(2), RBCOORDS(3)
1712:               ENDIF1706:               ENDIF
1713: 1707:               
1714:            ENDDO1708:            ENDDO
1715: 1709:            
1716:         ENDDO1710:         ENDDO
1717: 1711:         
1718:      ENDDO1712:      ENDDO
1719: 1713:      
1720:      RETURN1714:      RETURN
1721: 1715:      
1722:   ELSE IF (RIGID) THEN1716:   ELSE IF (RIGID) THEN
1723:      LUNIT=GETUNIT()1717:      LUNIT=GETUNIT()
1724:      OPEN(UNIT=LUNIT,FILE='rigid.xyz',STATUS='UNKNOWN')1718:      OPEN(UNIT=LUNIT,FILE='rigid.xyz',STATUS='UNKNOWN')
1725:      DO J1=1,NSAVE1719:      DO J1=1,NSAVE
1726:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES1720:         WRITE(LUNIT,'(I6)') (NATOMS/2)*NRBSITES
1727:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)1721:         WRITE(LUNIT,10) J1, QMIN(J1), FF(J1)
1728:         DO J2=1,NATOMS/21722:         DO J2=1,NATOMS/2
1729:            CALL RBIO(QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&1723:            CALL RBIO(QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3),&
1730:                 QMINP(J1,3*(NATOMS/2+J2-1)+1),&1724:                 QMINP(J1,3*(NATOMS/2+J2-1)+1),&
1731:                 QMINP(J1,3*(NATOMS/2+J2-1)+2),&1725:                 QMINP(J1,3*(NATOMS/2+J2-1)+2),&
1732:                 QMINP(J1,3*(NATOMS/2+J2-1)+3),&1726:                 QMINP(J1,3*(NATOMS/2+J2-1)+3),&
1733:                 RBCOORDS,NRBSITES,SITE)1727:                 RBCOORDS,NRBSITES,SITE)
1734:            DO J3=1,NRBSITES1728:            DO J3=1,NRBSITES
1735:               WRITE(LUNIT,'(A4,3F20.10)') 'LA ',RBCOORDS(3*(J3-1)+1),RBCOORDS(3*(J3-1)+2),RBCOORDS(3*(J3-1)+3)1729:               WRITE(LUNIT,'(A4,3F20.10)') 'LA ',RBCOORDS(3*(J3-1)+1),RBCOORDS(3*(J3-1)+2),RBCOORDS(3*(J3-1)+3)
1736:            ENDDO1730:            ENDDO
1737:         ENDDO1731:         ENDDO
1738:      ENDDO1732:      ENDDO
1739:      CLOSE(LUNIT)1733:      CLOSE(LUNIT)
1740: 1734:      
1741:   ENDIF1735:   ENDIF
1742: 1736:   
1743:   IF(ALLOCATED(DBNAME)) DEALLOCATE(DBNAME)1737:   IF(ALLOCATED(DBNAME)) DEALLOCATE(DBNAME)
1744: 1738:   
1745:   IF (AMBER12T) THEN1739:   IF (AMBER12T) THEN
1746:      CALL AMBER12_FINISH()1740:      CALL AMBER12_FINISH()
1747:   END IF1741:   END IF
1748: 1742:   
1749:   IF (OPEPT) CALL OPEP_FINISH()1743:   IF (OPEPT) CALL OPEP_FINISH()
1750: 1744:   
1751:   CALL CPU_TIME(TEND)1745:   CALL CPU_TIME(TEND)
1752:   WRITE(MYUNIT,"(A,F18.1,A)") "time elapsed ", TEND - TSTART, " seconds"1746:   WRITE(MYUNIT,"(A,F18.1,A)") "time elapsed ", TEND - TSTART, " seconds"
1753:   WRITE(MYUNIT,"(A,I18)") "Number of potential calls ", NPCALL1747:   WRITE(MYUNIT,"(A,I18)") "Number of potential calls ", NPCALL
1754:   RETURN1748:   RETURN
1755: 1749:   
1756: END SUBROUTINE FINALIO1750: END SUBROUTINE FINALIO
1757: 1751: 
1758: SUBROUTINE AMBERDUMP(J1,QMINP)1752: SUBROUTINE AMBERDUMP(J1,QMINP)
1759:     USE COMMONS1753:     USE COMMONS
1760:     USE MODAMBER1754:     USE MODAMBER
1761:     IMPLICIT NONE1755:     IMPLICIT NONE
1762: 1756: 
1763: 1757: 
1764:     CHARACTER(LEN=25) COORDFILE1758:     CHARACTER(LEN=25) COORDFILE
1765:     CHARACTER(LEN=2) FNAME1759:     CHARACTER(LEN=2) FNAME
1773:     ENDIF1767:     ENDIF
1774: 1768: 
1775:     DO A=1,ATOMS1769:     DO A=1,ATOMS
1776:         X(A)=QMINP(J1,3*A-2)1770:         X(A)=QMINP(J1,3*A-2)
1777:         Y(A)=QMINP(J1,3*A-1)1771:         Y(A)=QMINP(J1,3*A-1)
1778:         Z(A)=QMINP(J1,3*A)1772:         Z(A)=QMINP(J1,3*A)
1779:     END DO1773:     END DO
1780: 1774: 
1781:     COORDFILE='acoords.dump.'//FNAME1775:     COORDFILE='acoords.dump.'//FNAME
1782: 1776: 
1783: 1777:  
1784:     LUNIT=GETUNIT()1778:     LUNIT=GETUNIT()
1785:     OPEN (UNIT=LUNIT,IOSTAT=IOS,FILE=coordfile,STATUS='UNKNOWN')1779:     OPEN (UNIT=LUNIT,IOSTAT=IOS,FILE=coordfile,STATUS='UNKNOWN')
1786: 1780: 
1787:     DO a=1,atoms1781:     DO a=1,atoms
1788:         WRITE (UNIT=LUNIT,FMT='(A1,2X,A2,2X,I3,2X,I3,2X,F7.3,3X,F7.3,3X,F7.3)') label(a),typech(a),&1782:         WRITE (UNIT=LUNIT,FMT='(A1,2X,A2,2X,I3,2X,I3,2X,F7.3,3X,F7.3,3X,F7.3)') label(a),typech(a),&
1789:         a,bondedto(a),x(a),y(a),z(a)1783:         a,bondedto(a),x(a),y(a),z(a)
1790:     ENDDO1784:     ENDDO
1791: 1785: 
1792:     WRITE (UNIT=LUNIT,FMT='(A3)') 'end'1786:     WRITE (UNIT=LUNIT,FMT='(A3)') 'end'
1793:     WRITE (UNIT=LUNIT,FMT='(A)') ' '1787:     WRITE (UNIT=LUNIT,FMT='(A)') ' '
1882:     CA1=COS(ALPHA1)1876:     CA1=COS(ALPHA1)
1883:     C2A1=CA11877:     C2A1=CA1
1884:     IF (ALPHA1.LT.0.0001D0) THEN1878:     IF (ALPHA1.LT.0.0001D0) THEN
1885:         !        C3A1=(-ALPHA1/2+ALPHA1**3/24)1879:         !        C3A1=(-ALPHA1/2+ALPHA1**3/24)
1886:         C3A1=(-0.5D0+ALPHA1**2/24.0D0)1880:         C3A1=(-0.5D0+ALPHA1**2/24.0D0)
1887:         S1=(1.0D0-ALPHA1**2/6)1881:         S1=(1.0D0-ALPHA1**2/6)
1888:     ELSE1882:     ELSE
1889:         C3A1=(CA1-1.0D0)/ALPHA1**21883:         C3A1=(CA1-1.0D0)/ALPHA1**2
1890:         S1=SIN(ALPHA1)/ALPHA11884:         S1=SIN(ALPHA1)/ALPHA1
1891:     ENDIF1885:     ENDIF
1892: 1886:    
1893:     DO J1=1,NRBSITES1887:     DO J1=1,NRBSITES
1894:         COORDS(3*(J1-1)+1)=c2a1*SITE(J1,1) + s1*(n1*SITE(J1,2) - m1*SITE(J1,3)) - &1888:         COORDS(3*(J1-1)+1)=c2a1*SITE(J1,1) + s1*(n1*SITE(J1,2) - m1*SITE(J1,3)) - &
1895:         c3a1*l1*(l1*SITE(J1,1) + m1*SITE(J1,2) + n1*SITE(J1,3)) + X11889:         c3a1*l1*(l1*SITE(J1,1) + m1*SITE(J1,2) + n1*SITE(J1,3)) + X1
1896:         COORDS(3*(J1-1)+2)=c2a1*SITE(J1,2) + s1*(-(n1*SITE(J1,1)) + l1*SITE(J1,3)) &1890:         COORDS(3*(J1-1)+2)=c2a1*SITE(J1,2) + s1*(-(n1*SITE(J1,1)) + l1*SITE(J1,3)) &
1897:         - c3a1*m1*(l1*SITE(J1,1) + m1*SITE(J1,2) + n1*SITE(J1,3)) + Y11891:         - c3a1*m1*(l1*SITE(J1,1) + m1*SITE(J1,2) + n1*SITE(J1,3)) + Y1
1898:         COORDS(3*(J1-1)+3)=s1*(m1*SITE(J1,1) - l1*SITE(J1,2)) + c2a1*SITE(J1,3) &1892:         COORDS(3*(J1-1)+3)=s1*(m1*SITE(J1,1) - l1*SITE(J1,2)) + c2a1*SITE(J1,3) &
1899:         - c3a1*n1*(l1*SITE(J1,1) + m1*SITE(J1,2) + n1*SITE(J1,3)) + Z11893:         - c3a1*n1*(l1*SITE(J1,1) + m1*SITE(J1,2) + n1*SITE(J1,3)) + Z1
1900:     ENDDO1894:     ENDDO
1901: 1895: 
1902:     RETURN1896:     RETURN


r32062/keywords.f 2017-03-10 15:30:13.941255859 +0000 r32061/keywords.f 2017-03-10 15:30:14.729266279 +0000
 33:      &                      ligtransstep, ligmovefreq, amchnmax, amchnmin, amchpmax, amchpmin, rotamert, rotmaxchange, 33:      &                      ligtransstep, ligmovefreq, amchnmax, amchnmin, amchpmax, amchpmin, rotamert, rotmaxchange,
 34:      &                      rotcentre, rotpselect, rotoccuw, rotcutoff, setchiralgeneric, PRMTOP, IGB, RGBMAX, CUT, 34:      &                      rotcentre, rotpselect, rotoccuw, rotcutoff, setchiralgeneric, PRMTOP, IGB, RGBMAX, CUT,
 35:      &                      SALTCON, macroiont, nmacroions, macroiondist 35:      &                      SALTCON, macroiont, nmacroions, macroiondist
 36:       USE modamber 36:       USE modamber
 37:       USE PORFUNCS 37:       USE PORFUNCS
 38:       USE MYGA_PARAMS 38:       USE MYGA_PARAMS
 39:       USE BGUPMOD 39:       USE BGUPMOD
 40:       USE GLJYMOD 40:       USE GLJYMOD
 41:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L 41:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L
 42:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP 42:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP
 43:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, 43:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_RCUT, LJ_GAUSS_EPS, LJ_GAUSS_R0,
 44:      &                        LJ_GAUSS_R0, LJ_GAUSS_INITIALISE 44:      &                        LJ_GAUSS_INITIALISE
 45:       USE MBPOLMOD, ONLY: MBPOLINIT 45:       USE MBPOLMOD, ONLY: MBPOLINIT
 46:       USE SWMOD, ONLY: SWINIT, MWINIT 46:       USE SWMOD, ONLY: SWINIT, MWINIT
 47:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_GET_COORDS 47:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS,
 48: !     &                                 AMBER12_ATOMSS, AMBER12_SETUP, 48:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA
 49: !     &                                 AMBER12_RESIDUES, 
 50: !     &                                 POPULATE_ATOM_DATA 
 51:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 49:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 52:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 50:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 53:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 51:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,
 54:      &     PARSE_MLJ_PARAMS 52:      &     PARSE_MLJ_PARAMS
 55:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT 53:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT
 56:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE 54:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE
 57:       USE MOLECULAR_DYNAMICS, ONLY : MDT, MD_TSTEP, MD_GAMMA, MD_NWAIT, MD_NFREQ, MD_NSTEPS 55:       USE MOLECULAR_DYNAMICS, ONLY : MDT, MD_TSTEP, MD_GAMMA, MD_NWAIT, MD_NFREQ, MD_NSTEPS
 58:       USE OPEP_INTERFACE_MOD, ONLY : OPEP_INIT 56:       USE OPEP_INTERFACE_MOD, ONLY : OPEP_INIT
 59:  57: 
 60:       IMPLICIT NONE 58:       IMPLICIT NONE
 61:  59: 
 62:       DOUBLE PRECISION, ALLOCATABLE :: MLPMEAN(:), MLQMEAN(:) 60:       DOUBLE PRECISION, ALLOCATABLE :: MLPMEAN(:), MLQMEAN(:)
 63:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, LAST, IX, J1, JP, NPCOUNT, NDUMMY, J2, J3 61:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, LAST, IX, J1, JP, NPCOUNT, NDUMMY, INDEX, J2, J3, J4
 64:       INTEGER DATA_UNIT, FUNIT 62:       INTEGER DATA_UNIT, FUNIT
 65:       INTEGER MOVABLEATOMINDEX 63:       INTEGER MOVABLEATOMINDEX
 66:       LOGICAL CAT, YESNO, PERMFILE, CONFILE 64:       LOGICAL CAT, YESNO, PERMFILE, CONFILE
 67:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR, 65:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR,
 68:      &                NERROR, ECHO, LAST, CAT 66:      &                NERROR, ECHO, LAST, CAT
 69:        DOUBLE PRECISION XX, ROH, ROM, WTHETA 67:        DOUBLE PRECISION XX, ROH, ROM, WTHETA
 70:       LOGICAL END, SKIPBL, CLEAR, ECHO 68:       LOGICAL END, SKIPBL, CLEAR, ECHO
 71:       CHARACTER WORD*16,PBC*3,WORD2*10 69:       CHARACTER WORD*16,PBC*3,WORD2*10
 72:       DOUBLE PRECISION EAMLJA0, EAMLJBETA, EAMLJZ0, DUMMY 70:       DOUBLE PRECISION EAMLJA0, EAMLJBETA, EAMLJZ0, DUMMY
 73:       COMMON /EAMLJCOMM/ EAMLJA0, EAMLJBETA, EAMLJZ0 71:       COMMON /EAMLJCOMM/ EAMLJA0, EAMLJBETA, EAMLJZ0
 76:       COMMON /BSNEW/ SLENGTH, NOK, NBAD, EPS 74:       COMMON /BSNEW/ SLENGTH, NOK, NBAD, EPS
 77:       DOUBLE PRECISION EPS2, RAD, HEIGHT 75:       DOUBLE PRECISION EPS2, RAD, HEIGHT
 78:       COMMON /CAPS/ EPS2, RAD, HEIGHT 76:       COMMON /CAPS/ EPS2, RAD, HEIGHT
 79:  77: 
 80:       LOGICAL MLPDONE, MLPNORM, MLQDONE, MLQNORM 78:       LOGICAL MLPDONE, MLPNORM, MLQDONE, MLQNORM
 81: !     LOGICAL IGNOREBIN(HISTBINMAX), FIXBIN 79: !     LOGICAL IGNOREBIN(HISTBINMAX), FIXBIN
 82: !     COMMON /IG/ IGNOREBIN, FIXBIN 80: !     COMMON /IG/ IGNOREBIN, FIXBIN
 83: !      DOUBLE PRECISION    PMAX,PMIN,NMAX,NMIN,SIDESTEP 81: !      DOUBLE PRECISION    PMAX,PMIN,NMAX,NMIN,SIDESTEP
 84: !      COMMON /AMBWORD/    PMAX,PMIN,NMAX,NMIN,SIDESTEP 82: !      COMMON /AMBWORD/    PMAX,PMIN,NMAX,NMIN,SIDESTEP
 85:  83: 
 86:       INTEGER NATOM !, NDUM 84:       INTEGER NATOM, DMODE, NDUM, NONTYPEA
 87: ! 85: !
 88: ! These arrays should have dimension MXATMS 86: ! These arrays should have dimension MXATMS
 89: ! 87: !
 90:       DOUBLE PRECISION, ALLOCATABLE :: CHX(:), CHY(:), CHZ(:), CHMASS(:) 88:       DOUBLE PRECISION, ALLOCATABLE :: CHX(:), CHY(:), CHZ(:), CHMASS(:)
  89:       CHARACTER(LEN=1) DUMMYCH
 91:       CHARACTER(LEN=100) TOPFILE,PARFILE 90:       CHARACTER(LEN=100) TOPFILE,PARFILE
 92:       CHARACTER(LEN=20) UNSTRING 91:       CHARACTER(LEN=20) UNSTRING
 93:       DOUBLE PRECISION LJREPBB, LJATTBB, LJREPLL, LJATTLL, LJREPNN, LJATTNN, 92:       DOUBLE PRECISION LJREPBB, LJATTBB, LJREPLL, LJATTLL, LJREPNN, LJATTNN,
 94:      &                 HABLN, HBBLN, HCBLN, HDBLN, EABLN, EBBLN, ECBLN, EDBLN, TABLN, TBBLN, TCBLN, TDBLN 93:      &                 HABLN, HBBLN, HCBLN, HDBLN, EABLN, EBBLN, ECBLN, EDBLN, TABLN, TBBLN, TCBLN, TDBLN
 95:       DOUBLE PRECISION LJREPBL, LJATTBL, LJREPBN, LJATTBN, LJREPLN, LJATTLN 94:       DOUBLE PRECISION LJREPBL, LJATTBL, LJREPBN, LJATTBN, LJREPLN, LJATTLN
 96:  95: 
 97: !     DC430 > 96: !     DC430 >
 98:       DOUBLE PRECISION :: LPL, LPR 97:       DOUBLE PRECISION :: LPL, LPR
 99:       LOGICAL          :: RBSYMTEST     ! jdf43> 98:       LOGICAL          :: RBSYMTEST     ! jdf43>
100: ! 99: !
109: 108: 
110:       !ab2111> dihedral rotation stuff109:       !ab2111> dihedral rotation stuff
111:       INTEGER dihedraloffset,dihedralgroupsize,A1,A2,A3,A4110:       INTEGER dihedraloffset,dihedralgroupsize,A1,A2,A3,A4
112: 111: 
113: ! ab2111 > Reservoir stuff112: ! ab2111 > Reservoir stuff
114:       INTEGER VECLINE, VECNUM, RESCOUNT113:       INTEGER VECLINE, VECNUM, RESCOUNT
115: 114: 
116: ! hk286 - DAMPED GROUP MOVES115: ! hk286 - DAMPED GROUP MOVES
117:       DOUBLE PRECISION GROUPATT116:       DOUBLE PRECISION GROUPATT
118: 117: 
119: !      CHARACTER(LEN=120), DIMENSION(:,:), ALLOCATABLE :: KEY_WORDS118:       CHARACTER(LEN=120), DIMENSION(:,:), ALLOCATABLE :: KEY_WORDS
 119:       DOUBLE PRECISION DUMMY1(NATOMS)
120: 120: 
121:       INTEGER :: MAXNSETS121:       INTEGER :: MAXNSETS
122: 122: 
123:       CHARACTER(LEN=10) :: OPEP_DUMMY123:       CHARACTER(LEN=10) :: OPEP_DUMMY
124: 124: 
125: !      OPEN(5120, FILE = 'data')125: !      OPEN(5120, FILE = 'data')
126: !      CALL NEW_INPUT(5120, KEY_WORDS)126: !      CALL NEW_INPUT(5120, KEY_WORDS)
127: !      CLOSE(5120)127: !      CLOSE(5120)
128: !      PRINT *, KEY_WORDS128: !      PRINT *, KEY_WORDS
129: 129: 
1633:          STRESST=.TRUE.1633:          STRESST=.TRUE.
1634:          IF(NITEMS.GT.1) CALL READI(STRESS_MODE)1634:          IF(NITEMS.GT.1) CALL READI(STRESS_MODE)
1635: ! BLN1635: ! BLN
1636:       ELSE IF (WORD.EQ.'BLN') THEN1636:       ELSE IF (WORD.EQ.'BLN') THEN
1637:          BLNT=.TRUE.1637:          BLNT=.TRUE.
1638:          CALL READF(RK_R)1638:          CALL READF(RK_R)
1639:          CALL READF(RK_THETA)1639:          CALL READF(RK_THETA)
1640:          ALLOCATE(BEADLETTER(NATOMS),BLNSSTRUCT(NATOMS),1640:          ALLOCATE(BEADLETTER(NATOMS),BLNSSTRUCT(NATOMS),
1641:      &            LJREP_BLN(NATOMS,NATOMS),LJATT_BLN(NATOMS,NATOMS),A_BLN(NATOMS),B_BLN(NATOMS),C_BLN(NATOMS),D_BLN(NATOMS))1641:      &            LJREP_BLN(NATOMS,NATOMS),LJATT_BLN(NATOMS,NATOMS),A_BLN(NATOMS),B_BLN(NATOMS),C_BLN(NATOMS),D_BLN(NATOMS))
1642:          OPEN(UNIT=1,FILE='BLNsequence',STATUS='OLD')1642:          OPEN(UNIT=1,FILE='BLNsequence',STATUS='OLD')
1643:          READ(1,*) 1643:          READ(1,*) DUMMYCH
1644:          READ(1,*) LJREPBB, LJATTBB1644:          READ(1,*) LJREPBB, LJATTBB
1645:          READ(1,*) LJREPLL, LJATTLL1645:          READ(1,*) LJREPLL, LJATTLL
1646:          READ(1,*) LJREPNN, LJATTNN1646:          READ(1,*) LJREPNN, LJATTNN
1647:          READ(1,*)1647:          READ(1,*) DUMMYCH
1648:          READ(1,*)1648:          READ(1,*) DUMMYCH
1649:          READ(1,*) HABLN, HBBLN, HCBLN, HDBLN1649:          READ(1,*) HABLN, HBBLN, HCBLN, HDBLN
1650:          READ(1,*) EABLN, EBBLN, ECBLN, EDBLN1650:          READ(1,*) EABLN, EBBLN, ECBLN, EDBLN
1651:          READ(1,*) TABLN, TBBLN, TCBLN, TDBLN1651:          READ(1,*) TABLN, TBBLN, TCBLN, TDBLN
1652:          DO J1=1,NATOMS-11652:          DO J1=1,NATOMS-1
1653:             READ(1,'(A1)',ADVANCE='NO') BEADLETTER(J1)1653:             READ(1,'(A1)',ADVANCE='NO') BEADLETTER(J1)
1654:          ENDDO1654:          ENDDO
1655:          READ(1,'(A1)') BEADLETTER(NATOMS) ! this line is needed to advance the input line for the next read1655:          READ(1,'(A1)') BEADLETTER(NATOMS) ! this line is needed to advance the input line for the next read
1656:          DO J1=1,NATOMS-31656:          DO J1=1,NATOMS-3
1657:             READ(1,'(A1)',ADVANCE='NO') BLNSSTRUCT(J1)1657:             READ(1,'(A1)',ADVANCE='NO') BLNSSTRUCT(J1)
1658:          ENDDO1658:          ENDDO
1682:          CALL READF(RK_R)1682:          CALL READF(RK_R)
1683:          CALL READF(RK_THETA)1683:          CALL READF(RK_THETA)
1684:          IF (NITEMS.GT.3) THEN1684:          IF (NITEMS.GT.3) THEN
1685:             CALL READF(GOFACTOR)1685:             CALL READF(GOFACTOR)
1686:          ENDIF1686:          ENDIF
1687:          ALLOCATE(BEADLETTER(NATOMS),BLNSSTRUCT(NATOMS),1687:          ALLOCATE(BEADLETTER(NATOMS),BLNSSTRUCT(NATOMS),
1688:      &            LJREP_BLN(NATOMS,NATOMS),LJATT_BLN(NATOMS,NATOMS),A_BLN(NATOMS),B_BLN(NATOMS),C_BLN(NATOMS),D_BLN(NATOMS))1688:      &            LJREP_BLN(NATOMS,NATOMS),LJATT_BLN(NATOMS,NATOMS),A_BLN(NATOMS),B_BLN(NATOMS),C_BLN(NATOMS),D_BLN(NATOMS))
1689:          LJREP_BLN=01689:          LJREP_BLN=0
1690:          LJATT_BLN=01690:          LJATT_BLN=0
1691:          OPEN(UNIT=1,FILE='BLNsequence',STATUS='OLD')1691:          OPEN(UNIT=1,FILE='BLNsequence',STATUS='OLD')
1692:          READ(1,*)1692:          READ(1,*) DUMMYCH
1693:          READ(1,*) LJREPBB, LJATTBB1693:          READ(1,*) LJREPBB, LJATTBB
1694:          READ(1,*) LJREPLL, LJATTLL1694:          READ(1,*) LJREPLL, LJATTLL
1695:          READ(1,*) LJREPNN, LJATTNN1695:          READ(1,*) LJREPNN, LJATTNN
1696:          READ(1,*)1696:          READ(1,*) DUMMYCH
1697:          READ(1,*)1697:          READ(1,*) DUMMYCH
1698:          READ(1,*) HABLN, HBBLN, HCBLN, HDBLN1698:          READ(1,*) HABLN, HBBLN, HCBLN, HDBLN
1699:          READ(1,*) EABLN, EBBLN, ECBLN, EDBLN1699:          READ(1,*) EABLN, EBBLN, ECBLN, EDBLN
1700:          READ(1,*) TABLN, TBBLN, TCBLN, TDBLN1700:          READ(1,*) TABLN, TBBLN, TCBLN, TDBLN
1701:          DO J1=1,NATOMS-11701:          DO J1=1,NATOMS-1
1702:             READ(1,'(A1)',ADVANCE='NO') BEADLETTER(J1)1702:             READ(1,'(A1)',ADVANCE='NO') BEADLETTER(J1)
1703:          ENDDO1703:          ENDDO
1704:          READ(1,'(A1)') BEADLETTER(NATOMS) ! this line is needed to advance the input line for the next read1704:          READ(1,'(A1)') BEADLETTER(NATOMS) ! this line is needed to advance the input line for the next read
1705:          DO J1=1,NATOMS-31705:          DO J1=1,NATOMS-3
1706:             READ(1,'(A1)',ADVANCE='NO') BLNSSTRUCT(J1)1706:             READ(1,'(A1)',ADVANCE='NO') BLNSSTRUCT(J1)
1707:          ENDDO1707:          ENDDO
1729:          CHAPERONINT=.TRUE.1729:          CHAPERONINT=.TRUE.
1730:          CALL READF(RK_R)1730:          CALL READF(RK_R)
1731:          CALL READF(RK_THETA)1731:          CALL READF(RK_THETA)
1732:          CALL READF(RADIUS_CONTAINER)1732:          CALL READF(RADIUS_CONTAINER)
1733:          CALL READF(HYDROPHOBIC)1733:          CALL READF(HYDROPHOBIC)
1734:          ALLOCATE(BEADLETTER(NATOMS),BLNSSTRUCT(NATOMS),1734:          ALLOCATE(BEADLETTER(NATOMS),BLNSSTRUCT(NATOMS),
1735:      $        LJREP_BLN(NATOMS,NATOMS),LJATT_BLN(NATOMS,NATOMS)1735:      $        LJREP_BLN(NATOMS,NATOMS),LJATT_BLN(NATOMS,NATOMS)
1736:      $        ,A_BLN(NATOMS),B_BLN(NATOMS),C_BLN(NATOMS),D_BLN(NATOMS)1736:      $        ,A_BLN(NATOMS),B_BLN(NATOMS),C_BLN(NATOMS),D_BLN(NATOMS)
1737:      $        ,HYDRO_BLN(NATOMS))1737:      $        ,HYDRO_BLN(NATOMS))
1738:          OPEN(UNIT=1,FILE='BLNsequence',STATUS='OLD')1738:          OPEN(UNIT=1,FILE='BLNsequence',STATUS='OLD')
1739:          READ(1,*)1739:          READ(1,*) DUMMYCH
1740:          READ(1,*) LJREPBB, LJATTBB1740:          READ(1,*) LJREPBB, LJATTBB
1741:          READ(1,*) LJREPLL, LJATTLL1741:          READ(1,*) LJREPLL, LJATTLL
1742:          READ(1,*) LJREPNN, LJATTNN1742:          READ(1,*) LJREPNN, LJATTNN
1743:          READ(1,*) LJREPBL, LJATTBL1743:          READ(1,*) LJREPBL, LJATTBL
1744:          READ(1,*) LJREPBN, LJATTBN1744:          READ(1,*) LJREPBN, LJATTBN
1745:          READ(1,*) LJREPLN, LJATTLN1745:          READ(1,*) LJREPLN, LJATTLN
1746:          READ(1,*)1746:          READ(1,*) DUMMYCH
1747:          READ(1,*)1747:          READ(1,*) DUMMYCH
1748:          READ(1,*) HABLN, HBBLN, HCBLN, HDBLN1748:          READ(1,*) HABLN, HBBLN, HCBLN, HDBLN
1749:          READ(1,*) EABLN, EBBLN, ECBLN, EDBLN1749:          READ(1,*) EABLN, EBBLN, ECBLN, EDBLN
1750:          READ(1,*) TABLN, TBBLN, TCBLN, TDBLN1750:          READ(1,*) TABLN, TBBLN, TCBLN, TDBLN
1751:          DO J1=1,NATOMS-11751:          DO J1=1,NATOMS-1
1752:             READ(1,'(A1)',ADVANCE='NO') BEADLETTER(J1)1752:             READ(1,'(A1)',ADVANCE='NO') BEADLETTER(J1)
1753:          ENDDO1753:          ENDDO
1754:          READ(1,'(A1)') BEADLETTER(NATOMS) ! this line is needed to advance the input line for the next read1754:          READ(1,'(A1)') BEADLETTER(NATOMS) ! this line is needed to advance the input line for the next read
1755:          DO J1=1,NATOMS-31755:          DO J1=1,NATOMS-3
1756:             READ(1,'(A1)',ADVANCE='NO') BLNSSTRUCT(J1)1756:             READ(1,'(A1)',ADVANCE='NO') BLNSSTRUCT(J1)
1757:          ENDDO1757:          ENDDO
4441:                ALLOCATE(CONACTIVE(NCONSTRAINTFIX))4441:                ALLOCATE(CONACTIVE(NCONSTRAINTFIX))
4442:             ELSE4442:             ELSE
4443:                INQUIRE(FILE='congeom',EXIST=CONFILE)4443:                INQUIRE(FILE='congeom',EXIST=CONFILE)
4444:                NCONGEOM=04444:                NCONGEOM=0
4445:                IF (.NOT.CONFILE) THEN4445:                IF (.NOT.CONFILE) THEN
4446:                   WRITE(MYUNIT,'(A)') ' keyword> WARNING *** no congeom file found. Will use end point minima only.'4446:                   WRITE(MYUNIT,'(A)') ' keyword> WARNING *** no congeom file found. Will use end point minima only.'
4447:                ELSE4447:                ELSE
4448:                   LUNIT=GETUNIT()4448:                   LUNIT=GETUNIT()
4449:                   OPEN(LUNIT,FILE='congeom',STATUS='OLD')4449:                   OPEN(LUNIT,FILE='congeom',STATUS='OLD')
4450:                   DO4450:                   DO
4451:                      READ(LUNIT,*,END=864)4451:                      READ(LUNIT,*,END=864) DUMMY1(1)
4452:                      NCONGEOM=NCONGEOM+14452:                      NCONGEOM=NCONGEOM+1
4453:                   ENDDO4453:                   ENDDO
4454: 864               CONTINUE4454: 864               CONTINUE
4455:                   NCONGEOM=NCONGEOM/NATOMS4455:                   NCONGEOM=NCONGEOM/NATOMS
4456:                   ALLOCATE(CONGEOM(NCONGEOM,3*NATOMS))4456:                   ALLOCATE(CONGEOM(NCONGEOM,3*NATOMS))
4457:                   REWIND(LUNIT)4457:                   REWIND(LUNIT)
4458:                   DO J1=1,NCONGEOM4458:                   DO J1=1,NCONGEOM
4459:                      READ(LUNIT,*) CONGEOM(J1,1:3*NATOMS)4459:                      READ(LUNIT,*) CONGEOM(J1,1:3*NATOMS)
4460:                   ENDDO4460:                   ENDDO
4461:                   CLOSE(LUNIT)4461:                   CLOSE(LUNIT)
4486: ! constraints and get their geometries.4486: ! constraints and get their geometries.
4487: !4487: !
4488:             INQUIRE(FILE='congeom',EXIST=CONFILE)4488:             INQUIRE(FILE='congeom',EXIST=CONFILE)
4489:             NCONGEOM=04489:             NCONGEOM=0
4490:             IF (.NOT.CONFILE) THEN4490:             IF (.NOT.CONFILE) THEN
4491:                PRINT '(A)',' keyword> WARNING *** no congeom file found. Will use end point minima only.'4491:                PRINT '(A)',' keyword> WARNING *** no congeom file found. Will use end point minima only.'
4492:             ELSE4492:             ELSE
4493:                LUNIT=GETUNIT()4493:                LUNIT=GETUNIT()
4494:                OPEN(LUNIT,FILE='congeom',STATUS='OLD')4494:                OPEN(LUNIT,FILE='congeom',STATUS='OLD')
4495:                DO4495:                DO
4496:                   READ(LUNIT,*,END=863)4496:                   READ(LUNIT,*,END=863) DUMMY1(1)
4497:                   NCONGEOM=NCONGEOM+14497:                   NCONGEOM=NCONGEOM+1
4498:                ENDDO4498:                ENDDO
4499: 863            CONTINUE4499: 863            CONTINUE
4500:                NCONGEOM=NCONGEOM/NATOMS4500:                NCONGEOM=NCONGEOM/NATOMS
4501:                ALLOCATE(CONGEOM(NCONGEOM,3*NATOMS))4501:                ALLOCATE(CONGEOM(NCONGEOM,3*NATOMS))
4502:                REWIND(LUNIT)4502:                REWIND(LUNIT)
4503:                DO J1=1,NCONGEOM4503:                DO J1=1,NCONGEOM
4504:                   READ(LUNIT,*) CONGEOM(J1,1:3*NATOMS)4504:                   READ(LUNIT,*) CONGEOM(J1,1:3*NATOMS)
4505:                ENDDO4505:                ENDDO
4506:                CLOSE(LUNIT)4506:                CLOSE(LUNIT)
6190:          CALL READF(COULQ)6190:          CALL READF(COULQ)
6191:          CALL READF(COULSWAP)6191:          CALL READF(COULSWAP)
6192:          CALL READF(COULTEMP)6192:          CALL READF(COULTEMP)
6193:          NRBSITES=16193:          NRBSITES=1
6194:          ALLOCATE(SITE(NRBSITES,3))6194:          ALLOCATE(SITE(NRBSITES,3))
6195: !        Maybe the above two lines are not necessary!6195: !        Maybe the above two lines are not necessary!
6196: 6196: 
6197: ! jwrm2>6197: ! jwrm2>
6198:       ELSE IF (WORD.EQ.'LJ_GAUSS') THEN6198:       ELSE IF (WORD.EQ.'LJ_GAUSS') THEN
6199:           LJ_GAUSST = .TRUE.6199:           LJ_GAUSST = .TRUE.
6200:           CALL READI(LJ_GAUSS_MODE) 
6201:           CALL READF(LJ_GAUSS_RCUT)6200:           CALL READF(LJ_GAUSS_RCUT)
6202:           CALL READF(LJ_GAUSS_EPS)6201:           CALL READF(LJ_GAUSS_EPS)
6203:           CALL READF(LJ_GAUSS_R0)6202:           CALL READF(LJ_GAUSS_R0)
6204:           CALL LJ_GAUSS_INITIALISE()6203:           CALL LJ_GAUSS_INITIALISE()
6205: ! jwrm2> end6204: ! jwrm2> end
6206: 6205: 
6207:       ELSE IF (WORD.EQ.'STOCK') THEN6206:       ELSE IF (WORD.EQ.'STOCK') THEN
6208:          STOCKT=.TRUE.6207:          STOCKT=.TRUE.
6209:          RIGID=.TRUE.6208:          RIGID=.TRUE.
6210:          NRBSITES=16209:          NRBSITES=1


r32062/lj_gauss.f90 2017-03-10 15:30:14.205259350 +0000 r32061/lj_gauss.f90 2017-03-10 15:30:14.985269664 +0000
 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: ! 19: !
 20: ! This module provides a routine for calculating the energy and gradients with 20: ! This module provides a routine for calculating the energy and gradients with
 21: ! a double well Lennard-Jones + Gaussian potential. 21: ! a double well Lennard-Jones + Gaussian potential.
 22: ! Questions to jwrm2. 22: ! Questions to jwrm2.
 23: ! 23: !
 24: MODULE LJ_GAUSS_MOD 24: MODULE LJ_GAUSS_MOD
 25:  25: 
 26: PUBLIC :: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, LJ_GAUSS_R0 26: PUBLIC :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_RCUT
 27: PUBLIC :: LJ_GAUSS, LJ_GAUSS_INITIALISE, VIEW_LJ_GAUSS 27: PUBLIC :: LJ_GAUSS, LJ_GAUSS_INITIALISE
 28: PRIVATE :: LJ_GAUSS_SIGMASQ, RCUTINV6, RCUTINV12, EXP_RCUT, PREF_RCUT 28: PRIVATE :: LJ_GAUSS_SIGMASQ, RCUTINV6, RCUTINV12, EXP_RCUT, PREF_RCUT
 29: PRIVATE :: PERIODIC_RESET, LJ_GAUSS_PER, CALC_ENERGY, CALC_GRAD 29: PRIVATE :: PERIODIC_RESET, LJ_GAUSS_PER, CALC_ENERGY, CALC_GRAD
 30:  30: 
 31: ! EPS: strength of the Gaussian potential. The LJ strength is 1. 31: ! EPS: strength of the Gaussian potential. The LJ strength is 1.
 32: ! R0: Equilibrium distance of the minimum in the Gaussian potential. The LJ 32: ! R0: Equilibrium distance of the minimum in the Gaussian potential. The LJ
 33: !     rmin is 1. 33: !     rmin is 1.
 34: ! RCUTSQ: cutoff distance at which the energy, gradients and Hessian of 34: ! RCUTSQ: cutoff distance at which the energy, gradients and Hessian of
 35: !         both parts of the potential go smoothly to zero. 35: !         both parts of the potential go smoothly to zero.
 36: DOUBLE PRECISION :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_RCUT 36: DOUBLE PRECISION :: LJ_GAUSS_EPS, LJ_GAUSS_R0, LJ_GAUSS_RCUT
 37:  37: 
 38: ! LJ_GAUSS_MODE: 0 standard minimisation 
 39: !                1 the final triplet of coordinates will be interpreted as the 
 40: !                  unit cell parameters of an orthogonal unit cell, and 
 41: !                  minimised 
 42: INTEGER :: LJ_GAUSS_MODE 
 43:  
 44: ! LJ_GAUSS_SIGMA: variation (width) of the Gaussian potential 38: ! LJ_GAUSS_SIGMA: variation (width) of the Gaussian potential
 45: DOUBLE PRECISION, PARAMETER :: LJ_GAUSS_SIGMASQ = 0.02D0 39: DOUBLE PRECISION, PARAMETER :: LJ_GAUSS_SIGMASQ = 0.02D0
 46:  40: 
 47: ! RCUTINV6, RCUTINV12: factors for LJ which don't pened on particle distance 41: ! RCUTINV6, RCUTINV12: factors for LJ which don't pened on particle distance
 48: DOUBLE PRECISION :: RCUTINV6, RCUTINV12 42: DOUBLE PRECISION :: RCUTINV6, RCUTINV12
 49:  43: 
 50: ! EXP_RCUT, PREF_RCUT: factors for Gauss which don't depend on the particle 44: ! EXP_RCUT, PREF_RCUT: factors for Gauss which don't depend on the particle
 51: !                      distance 45: !                      distance
 52: DOUBLE PRECISION :: EXP_RCUT, PREF_RCUT 46: DOUBLE PRECISION :: EXP_RCUT, PREF_RCUT
 53:  47: 
 68:         ! PERIODIC: whether to use periodic boundary conditions 62:         ! PERIODIC: whether to use periodic boundary conditions
 69:         USE COMMONS, ONLY: NATOMS, PERIODIC 63:         USE COMMONS, ONLY: NATOMS, PERIODIC
 70:  64: 
 71:         IMPLICIT NONE 65:         IMPLICIT NONE
 72:  66: 
 73:         DOUBLE PRECISION, INTENT(INOUT) :: X(3*NATOMS) 67:         DOUBLE PRECISION, INTENT(INOUT) :: X(3*NATOMS)
 74:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), ENERGY 68:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOMS), ENERGY
 75:         LOGICAL, INTENT(IN) :: GTEST 69:         LOGICAL, INTENT(IN) :: GTEST
 76:  70: 
 77:         ! I, J: indices for looping over particles 71:         ! I, J: indices for looping over particles
 78:         ! NMOL: actual number of particles 72:         INTEGER :: I, J
 79:         INTEGER :: I, J, NMOL 
 80:  73: 
 81:         ! RIJ: Interparticle vector 74:         ! RIJ: Interparticle vector
 82:         ! MODRIJ: distance between a pair of particles 75:         ! MODRIJ: distance between a pair of particles
 83:         ! DVDR: derivative of pair potential wrt MODRIJ 76:         ! DVDR: derivative of pair potential wrt MODRIJ
 84:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR 77:         DOUBLE PRECISION :: RIJ(3), MODRIJ, DVDR
 85:  78: 
 86:         ! Initialise output variables 79:         ! Initialise output variables
 87:         ENERGY = 0.D0 80:         ENERGY = 0.D0
 88:         GRAD(:) = 0.D0 81:         GRAD(:) = 0.D0
 89:  82: 
 90:         IF (.NOT. PERIODIC .AND. LJ_GAUSS_MODE .EQ. 0) THEN 83:         ! Reset all particles to within the box
 91:             ! Normal cluster 84:         IF (PERIODIC) CALL PERIODIC_RESET(X(:))
 92:             NMOL = NATOMS 85: 
 93:             DO I = 1, NMOL-1 ! Outer loop over particles 86:         DO I = 1, NATOMS-1 ! Outer loop over particles
 94:                 DO J = I+1, NMOL ! Inner loop over particles 87:             DO J = I+1, NATOMS ! Inner loop over particles
  88: 
  89:                 IF (.NOT. PERIODIC) THEN
 95:                     ! Get the particle-particle distance 90:                     ! Get the particle-particle distance
 96:                     RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J) 91:                     RIJ(:) = X(3*I-2:3*I) - X(3*J-2:3*J)
 97:                     MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:))) 92:                     MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:)))
 98:  93: 
 99:                     ! Add energy and gradients 94:                     ! Add energy and gradients
100:                     IF (MODRIJ .LT. LJ_GAUSS_RCUT) THEN 95:                     IF (MODRIJ .LT. LJ_GAUSS_RCUT) THEN
101:                         ENERGY = ENERGY + CALC_ENERGY(MODRIJ) 96:                         ENERGY = ENERGY + CALC_ENERGY(MODRIJ)
102:  97: 
103:                         IF (GTEST) THEN 98:                         IF (GTEST) THEN
104:                             DVDR = CALC_GRAD(MODRIJ) 99:                             DVDR = CALC_GRAD(MODRIJ)
105:                             GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + &100:                             GRAD(3*I-2:3*I) = GRAD(3*I-2:3*I) + &
106:                                               DVDR*RIJ(:)/MODRIJ101:                                               DVDR*RIJ(:)/MODRIJ
107:                             GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - &102:                             GRAD(3*J-2:3*J) = GRAD(3*J-2:3*J) - &
108:                                               DVDR*RIJ(:)/MODRIJ103:                                               DVDR*RIJ(:)/MODRIJ
109:                         END IF ! GTEST104:                         END IF ! GTEST
110:                     END IF ! If less than cutoff distance105:                     END IF ! If less than cutoff distance
111:                 END DO ! Inner loop over particles 
112:             END DO ! Outer loop over particles 
113: 106: 
114:         ELSE IF (PERIODIC .AND. LJ_GAUSS_MODE .EQ. 0) THEN107:                 ELSE ! Periodic
115:             ! Periodic, fixed box 
116:             ! Reset all particles to within the box 
117:             CALL PERIODIC_RESET(X(:)) 
118:  
119:             NMOL = NATOMS 
120:             ! We need to include self-interactions of particles in different 
121:             ! unit cells 
122:             DO I = 1, NMOL ! Outer loop over particles 
123:                 DO J = I, NMOL ! Inner loop over particles 
124:                     ! Add energy and gradients108:                     ! Add energy and gradients
125:                     CALL LJ_GAUSS_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), &109:                     CALL LJ_GAUSS_PER(X(3*I-2:3*I), X(3*J-2:3*J), &
126:                                       GRAD(3*I-2:3*I), GRAD(3*J-2:3*J), &110:                                       GRAD(3*I-2:3*I), GRAD(3*J-2:3*J), &
127:                                       ENERGY, GTEST)111:                                       ENERGY, GTEST)
128:                 END DO ! Inner loop over particles112:                 END IF ! End if periodic
129:             END DO ! Outer loop over particles 
130: 113: 
131:         ELSE IF (PERIODIC .AND. LJ_GAUSS_MODE .EQ. 1) THEN114:             END DO ! Inner loop over particles
132:             ! Periodic varying lattice115:         END DO ! Outer loop over particles
133:             ! Reset all particles to within the box 
134:             CALL PERIODIC_RESET(X(:)) 
135:  
136:             NMOL = NATOMS-1 
137:             ! We need to include self-interactions of particles in different 
138:             ! unit cells 
139:             DO I = 1, NMOL! Outer loop over particles 
140:                 DO J = I, NMOL ! Inner loop over particles 
141:                     ! Add energy and gradients 
142:                     CALL LJ_GAUSS_PER(I, J, X(3*I-2:3*I), X(3*J-2:3*J), & 
143:                                       GRAD(3*I-2:3*I), GRAD(3*J-2:3*J), & 
144:                                       ENERGY, GTEST, GRAD(3*NATOMS-2:3*NATOMS)) 
145:                 END DO ! Inner loop over particles 
146:             END DO ! Outer loop over particles 
147:         END IF ! Mode selection 
148: 116: 
149:     END SUBROUTINE LJ_GAUSS117:     END SUBROUTINE LJ_GAUSS
150: 118: 
151: !-------------------------------------------------------------------------------119: !-------------------------------------------------------------------------------
152: !120: !
153: ! Calculates the energy contribution for a pair of particles.121: ! Calculates the energy contribution for a pair of particles.
154: ! MODRIJ: distance beteen the particle pair122: ! MODRIJ: distance beteen the particle pair
155: !123: !
156: !-------------------------------------------------------------------------------124: !-------------------------------------------------------------------------------
157:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ)125:     PURE DOUBLE PRECISION FUNCTION CALC_ENERGY(MODRIJ)
261: 229: 
262:     END FUNCTION CALC_GRAD230:     END FUNCTION CALC_GRAD
263: 231: 
264: !-------------------------------------------------------------------------------232: !-------------------------------------------------------------------------------
265: !233: !
266: ! Initialisation. Calculate some factors that won't be changing234: ! Initialisation. Calculate some factors that won't be changing
267: !235: !
268: !-------------------------------------------------------------------------------236: !-------------------------------------------------------------------------------
269:     SUBROUTINE LJ_GAUSS_INITIALISE ()237:     SUBROUTINE LJ_GAUSS_INITIALISE ()
270: 238: 
271:         ! MYUNIT: file unit for main GMIN output 
272:         USE COMMONS, ONLY: MYUNIT 
273:  
274:         IMPLICIT NONE239:         IMPLICIT NONE
275: 240: 
276:         ! Check we have sensible value for the mode 
277:         IF (LJ_GAUSS_MODE .NE. 1 .AND. LJ_GAUSS_MODE .NE. 0) THEN 
278:             WRITE (MYUNIT, *) 'LJ_GAUSS> ERROR: mode must be 0 or 1' 
279:             STOP 
280:         END IF 
281:  
282:         ! Calculate some factors for LJ that don't depend on RIJ241:         ! Calculate some factors for LJ that don't depend on RIJ
283:         RCUTINV6 = (1.D0/LJ_GAUSS_RCUT)**6242:         RCUTINV6 = (1.D0/LJ_GAUSS_RCUT)**6
284:         RCUTINV12 = RCUTINV6**2243:         RCUTINV12 = RCUTINV6**2
285: 244: 
286:         ! Calculate some factors for Gauss that don't depend on RIJ245:         ! Calculate some factors for Gauss that don't depend on RIJ
287:         EXP_RCUT = LJ_GAUSS_EPS * EXP(- (LJ_GAUSS_RCUT - LJ_GAUSS_R0)**2 / &246:         EXP_RCUT = LJ_GAUSS_EPS * EXP(- (LJ_GAUSS_RCUT - LJ_GAUSS_R0)**2 / &
288:                                         (2.D0*LJ_GAUSS_SIGMASQ))247:                                         (2.D0*LJ_GAUSS_SIGMASQ))
289:         PREF_RCUT = (LJ_GAUSS_RCUT - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ248:         PREF_RCUT = (LJ_GAUSS_RCUT - LJ_GAUSS_R0) / LJ_GAUSS_SIGMASQ
290: 249: 
291:     END SUBROUTINE LJ_GAUSS_INITIALISE250:     END SUBROUTINE LJ_GAUSS_INITIALISE
292: 251: 
293: !-------------------------------------------------------------------------------252: !-------------------------------------------------------------------------------
294: !253: !
295: ! Resets all coordinates to within the unit cell254: ! Resets all coordinates to within the unit cell
296: ! The unit cell runs from 0 to 1 in all directions. Coordinates are later scaled 
297: ! by the box lengths. 
298: ! COORDS: coordinates of the particles255: ! COORDS: coordinates of the particles
299: !256: !
300: !-------------------------------------------------------------------------------257: !-------------------------------------------------------------------------------
301:     SUBROUTINE PERIODIC_RESET(COORDS)258:     SUBROUTINE PERIODIC_RESET(COORDS)
302: 259: 
303:         ! BOXLX, BOXLY, BOXLZ: dimensions of box260:         ! BOXLX, BOXLY, BOXLZ: dimensions of box
 261:         ! MYUNIT: file unit for main GMIN output
304:         ! NATOMS: number of particles262:         ! NATOMS: number of particles
305:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ, NATOMS263:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ, MYUNIT, NATOMS
306: 264: 
307:         IMPLICIT NONE265:         IMPLICIT NONE
308: 266: 
309:         DOUBLE PRECISION, INTENT(INOUT) :: COORDS(3*NATOMS)267:         DOUBLE PRECISION, INTENT(INOUT) :: COORDS(3*NATOMS)
310: 268: 
311:         ! Sort out the unit cell sizes, if they are varying269:         INTEGER :: I
312:         IF (LJ_GAUSS_MODE .EQ. 1) THEN 
313:             BOXLX = COORDS(3*NATOMS - 2) 
314:             BOXLY = COORDS(3*NATOMS - 1) 
315:             BOXLZ = COORDS(3*NATOMS) 
316:         END IF 
317: 270: 
318:         ! Reset coordinates to within the box271:         ! Check RCUT against the box sizes
319:         IF (LJ_GAUSS_MODE .EQ. 0) THEN272:         IF (LJ_GAUSS_RCUT .GT. BOXLX .OR. LJ_GAUSS_RCUT .GT. BOXLY &
320:             COORDS(:) = COORDS(:) - FLOOR(COORDS(:))273:             .OR. LJ_GAUSS_RCUT .GT. BOXLZ) THEN
321:         ELSE IF (LJ_GAUSS_MODE .EQ. 1) THEN274:             WRITE (MYUNIT, *) 'LJ_GAUSS> Box size less the potential cutoff.', &
322:             COORDS(1:3*NATOMS-3) = COORDS(1:3*NATOMS-3) - &275:                               ' This can cause discontinuities.'
323:                                    FLOOR(COORDS(1:3*NATOMS-3)) 
324:         END IF276:         END IF
325: 277: 
 278:         DO I = 1, NATOMS ! Loop over particles
 279:             ! The box stretches from -0.5*BOXLx to 0.5*BOXLx in each direction
 280:             COORDS(3*I-2) = COORDS(3*I - 2) - BOXLX*NINT(COORDS(3*I - 2)/BOXLX)
 281:             COORDS(3*I-1) = COORDS(3*I - 1) - BOXLY*NINT(COORDS(3*I - 1)/BOXLY)
 282:             COORDS(3*I  ) = COORDS(3*I    ) - BOXLZ*NINT(COORDS(3*I    )/BOXLZ)
 283:         END DO ! Loop over particles
 284: 
326:     END SUBROUTINE PERIODIC_RESET285:     END SUBROUTINE PERIODIC_RESET
327: 286: 
328: !-------------------------------------------------------------------------------287: !-------------------------------------------------------------------------------
329: !288: !
330: ! Finds the energy and gradients for the periodic system289: ! Finds the energy and gradients for the periodic system
331: ! ID1, ID2: id numbers of the particles 
332: ! COORDS_I: coordinates of particle I290: ! COORDS_I: coordinates of particle I
333: ! COORDS_J: coordinates of particle J291: ! COORDS_J: coordinates of particle J
334: ! GRAD_I: gradients for particle I292: ! GRAD_I: gradients for particle I
335: ! GRAD_J: gradients for particle J293: ! GRAD_J: gradients for particle J
336: ! ENERGY: energy contribution for this pair294: ! ENERGY: energy contribution for this pair
337: ! RIJ: mininum length vector between the particles295: ! RIJ: mininum length vector between the particles
338: ! GTEST: whether to calculate gradients296: ! GTEST: whether to calculate gradients
339: ! GRAD_BOX: gradients for the box 
340: !297: !
341: !-------------------------------------------------------------------------------298: !-------------------------------------------------------------------------------
342: 299: 
343:     SUBROUTINE LJ_GAUSS_PER(ID1, ID2, COORDS_I,COORDS_J,GRAD_I,GRAD_J,ENERGY, &300:     SUBROUTINE LJ_GAUSS_PER(COORDS_I,COORDS_J,GRAD_I,GRAD_J,ENERGY,GTEST)
344:                             GTEST, GRAD_BOX) 
345: 301: 
346:         ! BOXLX, BOXLY, BOXLZ: dimensions of the box302:         ! BOXLX, BOXLY, BOXLZ: dimensions of the box
347:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ303:         USE COMMONS, ONLY: BOXLX, BOXLY, BOXLZ
348: 304: 
349:         INTEGER, INTENT(IN) :: ID1, ID2 
350:         DOUBLE PRECISION, INTENT(IN) :: COORDS_I(3), COORDS_J(3)305:         DOUBLE PRECISION, INTENT(IN) :: COORDS_I(3), COORDS_J(3)
351:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY306:         DOUBLE PRECISION, INTENT(INOUT) :: GRAD_I(3), GRAD_J(3), ENERGY
352:         DOUBLE PRECISION, INTENT(INOUT), OPTIONAL :: GRAD_BOX(3) 
353:         LOGICAL, INTENT(IN) :: GTEST307:         LOGICAL, INTENT(IN) :: GTEST
354: 308: 
355:         ! XJ: copy of j coordinates309:         ! XJ: copy of j coordinates
356:         ! RIJ(3): particle-particle vector310:         ! RIJ(3): particle-particle vector
357:         ! MODRIJ: particle-particle distance311:         ! MODRIJ: particle-particle distance
358:         ! DVDR: derivative of pair potential wrt MODRIJ312:         ! DVDR: derivative of pair potential wrt MODRIJ
359:         ! BOX_SIZE: convenience copy of the box size313:         DOUBLE PRECISION :: XJ(3), RIJ(3), MODRIJ, DVDR
360:         DOUBLE PRECISION :: XJ(3), RIJ(3), MODRIJ, DVDR, BOX_SIZE(3) 
361: 314: 
362:         ! PER_k: integers for looping over cell possibilities315:         ! PER: integer for looping over cell possibilities
363:         ! CELL_RANGE: number of cells in each direction required316:         INTEGER :: PER
364:         INTEGER :: PER_X, PER_Y, PER_Z, CELL_RANGE(3) 
365:  
366: !        DOUBLE PRECISION :: START_ENERGY 
367: !        START_ENERGY = ENERGY 
368:  
369:         BOX_SIZE(1) = BOXLX 
370:         BOX_SIZE(2) = BOXLY 
371:         BOX_SIZE(3) = BOXLZ 
372:  
373:         ! Set the number of cells we need to check, based on the cell size and 
374:         ! the potential cutoff. 
375:         CELL_RANGE(:) = CEILING(LJ_GAUSS_RCUT/BOX_SIZE(:)) 
376: 317: 
377:         ! Loop over different coordinate possibilities.318:         ! Loop over different coordinate possibilities.
378:         DO PER_X = - CELL_RANGE(1), CELL_RANGE(1)319:         DO PER = 0, 26
379:             DO PER_Y = - CELL_RANGE(2), CELL_RANGE(2) 
380:                 DO PER_Z = - CELL_RANGE(3), CELL_RANGE(3) 
381:  
382:                     ! Skip the interaction of a particle with itself in the same 
383:                     ! cell and don't double count interactions within the unit 
384:                     ! cell. 
385:                     IF (PER_X .EQ. 0 .AND. PER_Y .EQ. 0 .AND. PER_Z .EQ. 0 & 
386:                         .AND. ID1 .EQ. ID2) CYCLE 
387:  
388:                     ! Adjust the j coordinates for the periodic image 
389:                     XJ(1) = COORDS_J(1) + PER_X 
390:                     XJ(2) = COORDS_J(2) + PER_Y 
391:                     XJ(3) = COORDS_J(3) + PER_Z 
392: 320: 
393:                     ! Find the distance and scale by the box size321:             XJ(:) = COORDS_J(:)
394:                     RIJ(:) = (COORDS_I(:) - XJ(:))*BOX_SIZE(:) 
395:                     MODRIJ = SQRT(DOT_PRODUCT(RIJ(:), RIJ(:))) 
396: 322: 
397:                     IF (MODRIJ .LT. LJ_GAUSS_RCUT) THEN323:             XJ(1) = XJ(1) + (MOD(PER, 3) - 1) * BOXLX
398:                         ! Add energy and gradients324:             XJ(2) = XJ(2) + (MOD(PER/3, 3) - 1) * BOXLY
399:                         ! Divide by 2 for images of the same atom, to avoid325:             XJ(3) = XJ(3) + (MOD(PER/9, 3) - 1) * BOXLZ
400:                         ! double counting326: 
401:                         ENERGY = ENERGY + MERGE(0.5D0, 1.0D0, ID1 .EQ. ID2) * &327:             ! Find the distance.
402:                                           CALC_ENERGY(MODRIJ)328:             RIJ(:) = COORDS_I(:) - XJ(:)
 329:             MODRIJ = SQRT(DOT_PRODUCT(RIJ(:),RIJ(:)))
 330: 
 331:             IF (MODRIJ .LT. LJ_GAUSS_RCUT) THEN
 332:                 ! Add energy and gradients
 333:                 ENERGY = ENERGY + CALC_ENERGY(MODRIJ)
 334: 
 335:                 IF (GTEST) THEN
 336:                     DVDR = CALC_GRAD(MODRIJ)
 337:                     GRAD_I(:) = GRAD_I(:) + DVDR*RIJ(:)/MODRIJ
 338:                     GRAD_J(:) = GRAD_J(:) - DVDR*RIJ(:)/MODRIJ
 339:                 END IF ! GTEST
 340:             END IF ! End if less than cutoff
403: 341: 
404:                         IF (GTEST) THEN342:         END DO
405:                             ! Divide gradient by 2 for images of the same atom, 
406:                             ! to avoid double counting 
407:                             DVDR = MERGE(0.5D0, 1.0D0, ID1 .EQ. ID2) * & 
408:                                    CALC_GRAD(MODRIJ) 
409:                             ! We must undo the box size coordinate scaling 
410:                             GRAD_I(:) = GRAD_I(:) + & 
411:                                         DVDR*RIJ(:)*BOX_SIZE(:)/MODRIJ 
412:                             GRAD_J(:) = GRAD_J(:) - & 
413:                                         DVDR*RIJ(:)*BOX_SIZE(:)/MODRIJ 
414:  
415:                             IF (PRESENT(GRAD_BOX)) THEN 
416:                                 ! Lattice derivatives 
417:                                 GRAD_BOX(:) = GRAD_BOX(:) + & 
418:                                          DVDR*RIJ(:)*RIJ(:)/(BOX_SIZE(:)*MODRIJ) 
419:                             END IF ! lattice derivatievs 
420:                         END IF ! GTEST 
421:                     END IF ! End if less than cutoff 
422:                 END DO ! z loop 
423:             END DO ! y loop 
424:         END DO ! x loop 
425: 343: 
426:     END SUBROUTINE LJ_GAUSS_PER344:     END SUBROUTINE LJ_GAUSS_PER
427: 345: 
428: !------------------------------------------------------------------------------- 
429: ! 
430: ! Outputs the coordinates to the file ljgauss.xyz 
431: ! 
432: !------------------------------------------------------------------------------- 
433:     SUBROUTINE VIEW_LJ_GAUSS() 
434:  
435:         ! NATOMS: number of coordinates 
436:         ! NSAVE: number of minima to write 
437:         ! BOXx: Box dimensions 
438:         ! PERIODIC: whether periodic boundary conditions are in use 
439:         USE COMMONS, ONLY: NATOMS, NSAVE, PERIODIC, BOXLX, BOXLY, BOXLZ 
440:  
441:         ! QMIN: energies of the saved minima 
442:         ! QMINP: coordinates of the saved minima 
443:         ! FF: step at which the saved minima was first found 
444:         USE QMODULE, ONLY: QMIN, QMINP, FF 
445:  
446:         IMPLICIT NONE 
447:  
448:         ! GETUNIT: function for fiding a free file unit 
449:         ! FUNIT: output file unit 
450:         ! NMOL: actual number of particles 
451:         ! I, J: loop iterators 
452:         INTEGER :: GETUNIT, FUNIT, NMOL, I, J 
453:  
454:         ! BOX_SIZE: convenience array for the box sizes 
455:         DOUBLE PRECISION :: BOX_SIZE(3) 
456:  
457:         ! Set NMOL and box parameters 
458:         IF (LJ_GAUSS_MODE .EQ. 1) THEN 
459:             NMOL = NATOMS - 1 
460:         ELSE IF (LJ_GAUSS_MODE .EQ. 0 .AND. PERIODIC) THEN 
461:             NMOL = NATOMS 
462:             BOX_SIZE(1) = BOXLX 
463:             BOX_SIZE(2) = BOXLY 
464:             BOX_SIZE(3) = BOXLZ 
465:         ELSE 
466:             NMOL = NATOMS 
467:         END IF 
468:  
469:         ! Open the file 
470:         FUNIT = GETUNIT() 
471:         OPEN(UNIT=FUNIT, FILE='ljgauss.xyz', STATUS='UNKNOWN') 
472:  
473:         ! Loop over the configurations 
474:         DO I = 1, NSAVE 
475:  
476:             WRITE(FUNIT,'(I6)') NMOL 
477:  
478:             WRITE(FUNIT, '(A,I6,A,F20.10,A,I8,A,F20.10)') 'Energy of minimum ',& 
479:                 I, '=', QMIN(I), ' first found at step ', FF(I), & 
480:                 ' Energy per particle = ', QMIN(I)/NMOL 
481:  
482:             ! Set box parameters for varying box 
483:             IF (LJ_GAUSS_MODE .EQ. 1) THEN 
484:                 BOX_SIZE(:) = QMINP(I, 3*NATOMS-2:3*NATOMS) 
485:             END IF 
486:  
487:             DO J = 1, NMOL 
488:  
489:                 IF (PERIODIC) THEN 
490:                     IF (J .EQ. 1) THEN 
491:                         WRITE(FUNIT, '(A,3F20.10,A,3(A,F20.10))') 'O ', & 
492:                               QMINP(I, 3*J-2:3*J)*BOX_SIZE(:), ' bbox_xyz', & 
493:                               ' 0.0', BOX_SIZE(1), ' 0.0', BOX_SIZE(2), & 
494:                               ' 0.0', BOX_SIZE(3) 
495:                     ELSE 
496:                         WRITE(FUNIT, '(A,3F20.10)') 'O ', & 
497:                             QMINP(I, 3*J-2:3*J)*BOX_SIZE(:) 
498:                     END IF 
499:                 ELSE 
500:                     WRITE(FUNIT, '(A,3F20.10)') 'O ', QMINP(I, 3*J-2:3*J) 
501:                 END IF 
502:  
503:             END DO ! Loop over particles 
504:  
505:         END DO ! Loop over the configurations 
506:  
507:         CLOSE(FUNIT) 
508:  
509:     END SUBROUTINE VIEW_LJ_GAUSS 
510:  
511: END MODULE LJ_GAUSS_MOD346: END MODULE LJ_GAUSS_MOD


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0