hdiff output

r31387/commons.f90 2016-10-26 09:30:14.220627916 +0100 r31386/commons.f90 2016-10-26 09:30:16.168654020 +0100
319:       INTEGER :: HOMOREF_FGMODE, HOMOREF_LSMODE, HOMOREF_NCYCLES319:       INTEGER :: HOMOREF_FGMODE, HOMOREF_LSMODE, HOMOREF_NCYCLES
320:       INTEGER :: HOMOREF_AUX_NSWAPS, HOMOREF_NFMAX, HOMOREF_NSMAX320:       INTEGER :: HOMOREF_AUX_NSWAPS, HOMOREF_NFMAX, HOMOREF_NSMAX
321:       INTEGER :: HOMOREF_BH_NSWAPMAX, HOMOREF_BH_NDUDMAX321:       INTEGER :: HOMOREF_BH_NSWAPMAX, HOMOREF_BH_NDUDMAX
322:       INTEGER :: RANDPERM_STEP, RANDMULTIPERM_STEP322:       INTEGER :: RANDPERM_STEP, RANDMULTIPERM_STEP
323:       DOUBLE PRECISION :: HOMOREF_AUX_TEMP, HOMOREF_AUX_FACTOR323:       DOUBLE PRECISION :: HOMOREF_AUX_TEMP, HOMOREF_AUX_FACTOR
324:       DOUBLE PRECISION :: HOMOREF_AUX_NNCUT324:       DOUBLE PRECISION :: HOMOREF_AUX_NNCUT
325:       DOUBLE PRECISION :: HOMOREF_BH_TEMP, HOMOREF_BH_FACTOR325:       DOUBLE PRECISION :: HOMOREF_BH_TEMP, HOMOREF_BH_FACTOR
326: !ds656> Quench-Assisted Local Combinatorial Search326: !ds656> Quench-Assisted Local Combinatorial Search
327:       LOGICAL :: QALCST, QALCSV, QALCS_SURFT, QALCS_SYMT, SPECLABELST327:       LOGICAL :: QALCST, QALCSV, QALCS_SURFT, QALCS_SYMT, SPECLABELST
328:       INTEGER :: QALCSMODE, QALCS_SYM_MINCORESIZE, QALCS_NBRHD, &328:       INTEGER :: QALCSMODE, QALCS_SYM_MINCORESIZE, QALCS_NBRHD, &
329:            QALCS_PARAM, SEQLENGTH, QALCS_SURF_MODE329:            QALCS_PARAM, SEQLENGTH
330:       CHARACTER(LEN=2), ALLOCATABLE :: SPECLABELS(:)330:       CHARACTER(LEN=2), ALLOCATABLE :: SPECLABELS(:)
331: !ds656> Enumerate permutations331: !ds656> Enumerate permutations
332:       LOGICAL :: ENPERMST, MULTIPERMT, SPANSWAPST332:       LOGICAL :: ENPERMST, MULTIPERMT, SPANSWAPST
333:       DOUBLE PRECISION, ALLOCATABLE :: LEHMER_COORDS(:,:) ! NAT, NPAR333:       DOUBLE PRECISION, ALLOCATABLE :: LEHMER_COORDS(:,:) ! NAT, NPAR
334:       CHARACTER(LEN=1), ALLOCATABLE :: LEHMER_LIST(:,:) ! NATOMS, NPAR334:       CHARACTER(LEN=1), ALLOCATABLE :: LEHMER_LIST(:,:) ! NATOMS, NPAR
335:       CHARACTER(LEN=1), ALLOCATABLE :: LEHMER_LAST(:) ! NPAR335:       CHARACTER(LEN=1), ALLOCATABLE :: LEHMER_LAST(:) ! NPAR
336:       INTEGER, ALLOCATABLE :: LEHMER_ILASTB(:) ! NPAR336:       INTEGER, ALLOCATABLE :: LEHMER_ILASTB(:) ! NPAR
337: !ds656> Generalized LJ + Yukawa potential337: !ds656> Generalized LJ + Yukawa potential
338:       LOGICAL :: GLJY338:       LOGICAL :: GLJY
339: 339: 


r31387/keywords.f 2016-10-26 09:30:14.472631293 +0100 r31386/keywords.f 2016-10-26 09:30:16.432657563 +0100
496:       LFLIPS_TFACTOR = 1.0496:       LFLIPS_TFACTOR = 1.0
497:       LFLIPS_NUP = 1497:       LFLIPS_NUP = 1
498:       NTYPEA = NATOMS ! ds656> before this was initialised to 0498:       NTYPEA = NATOMS ! ds656> before this was initialised to 0
499:       BGUPTANAME1 = 'Au'499:       BGUPTANAME1 = 'Au'
500:       BGUPTANAME2 = 'Ag'500:       BGUPTANAME2 = 'Ag'
501:       QALCST = .FALSE.501:       QALCST = .FALSE.
502:       QALCSV = .FALSE.502:       QALCSV = .FALSE.
503:       QALCSMODE = 0503:       QALCSMODE = 0
504:       QALCS_PARAM = 0504:       QALCS_PARAM = 0
505:       QALCS_SURFT = .FALSE.505:       QALCS_SURFT = .FALSE.
506:       QALCS_SURF_MODE = 0 
507:       QALCS_SYMT = .FALSE.506:       QALCS_SYMT = .FALSE.
508:       QALCS_SYM_MINCORESIZE = 5507:       QALCS_SYM_MINCORESIZE = 5
509:       SPECLABELST = .FALSE.508:       SPECLABELST = .FALSE.
510:       RANDMULTIPERMT = .FALSE.509:       RANDMULTIPERMT = .FALSE.
511:       RANDMULTIPERM_STEP=1510:       RANDMULTIPERM_STEP=1
512:       RANDPERMT = .FALSE.511:       RANDPERMT = .FALSE.
513:       RANDPERM_STEP=1512:       RANDPERM_STEP=1
514:       MULTIPERMT = .FALSE.513:       MULTIPERMT = .FALSE.
515:       SPANSWAPST = .FALSE.514:       SPANSWAPST = .FALSE.
516:       ENPERMST = .FALSE.515:       ENPERMST = .FALSE.
1471:             CALL READI(QALCSMODE)1470:             CALL READI(QALCSMODE)
1472:             IF(NITEMS > 2) THEN1471:             IF(NITEMS > 2) THEN
1473:                CALL READI(QALCS_PARAM)1472:                CALL READI(QALCS_PARAM)
1474:             ENDIF1473:             ENDIF
1475:          ELSE1474:          ELSE
1476:             WRITE(MYUNIT,'(A)') 'keywords> QALCS mode unspecified!'1475:             WRITE(MYUNIT,'(A)') 'keywords> QALCS mode unspecified!'
1477:             STOP1476:             STOP
1478:          ENDIF1477:          ENDIF
1479:       ELSE IF (WORD .EQ. 'QALCS_SURF') THEN1478:       ELSE IF (WORD .EQ. 'QALCS_SURF') THEN
1480:          QALCS_SURFT = .TRUE.1479:          QALCS_SURFT = .TRUE.
1481:          IF(NITEMS > 1) CALL READI(QALCS_SURF_MODE) 
1482:       ELSE IF (WORD .EQ. 'QALCS_SYM') THEN1480:       ELSE IF (WORD .EQ. 'QALCS_SYM') THEN
1483:          QALCS_SYMT = .TRUE.1481:          QALCS_SYMT = .TRUE.
1484:          IF(NITEMS > 1) THEN1482:          IF(NITEMS > 1) THEN
1485:             CALL READI(QALCS_SYM_MINCORESIZE)1483:             CALL READI(QALCS_SYM_MINCORESIZE)
1486:          ENDIF1484:          ENDIF
1487:       ELSE IF (WORD .EQ. 'QALCSVERBOSE') THEN1485:       ELSE IF (WORD .EQ. 'QALCSVERBOSE') THEN
1488:          QALCSV = .TRUE.1486:          QALCSV = .TRUE.
1489:       ELSE IF(WORD .EQ. 'SPECLABELS') THEN1487:       ELSE IF(WORD .EQ. 'SPECLABELS') THEN
1490:          SPECLABELST = .TRUE.1488:          SPECLABELST = .TRUE.
1491:          ALLOCATE(SPECLABELS(NITEMS-1))1489:          ALLOCATE(SPECLABELS(NITEMS-1))


r31387/mc.F 2016-10-26 09:30:14.724634670 +0100 r31386/mc.F 2016-10-26 09:30:16.684660932 +0100
1547:                IF ( LSWAPST .AND.1547:                IF ( LSWAPST .AND.
1548:      1              (MOD(J1,LSWAPS_N1)==0) ) THEN1548:      1              (MOD(J1,LSWAPS_N1)==0) ) THEN
1549:                   CALL USWAPS(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)1549:                   CALL USWAPS(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
1550:                ENDIF1550:                ENDIF
1551: 1551: 
1552:                IF ( LFLIPST .AND.1552:                IF ( LFLIPST .AND.
1553:      1              (MOD(J1,LFLIPS_N1)==0) ) THEN1553:      1              (MOD(J1,LFLIPS_N1)==0) ) THEN
1554:                   CALL FLIPSEQ(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)1554:                   CALL FLIPSEQ(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
1555:                ENDIF1555:                ENDIF
1556: 1556: 
1557:                IF ( (QDONE.EQ.1) .AND. (.NOT.SYMMETRIZE) .AND.1557:                IF ((QDONE.EQ.1).AND.(.NOT.SYMMETRIZE)) THEN
1558:      1              (ABS(POTEL - EPREV(JP)) > ECONV .OR. NQ(JP)==1) ) THEN 
1559:                   IF (QALCST.OR.QALCS_SURFT.OR.QALCS_SYMT) CALL QALCS(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)1558:                   IF (QALCST.OR.QALCS_SURFT.OR.QALCS_SYMT) CALL QALCS(JP,ITERATIONS,TIME,BRUN,QDONE,SCREENC)
1560:                ENDIF1559:                ENDIF
1561: !<ds6561560: !<ds656
1562:                1561:                
1563: !  RMS compared to reference structure 'compare'1562: !  RMS compared to reference structure 'compare'
1564:           IF (RMST.AND.CHRMMT) THEN1563:           IF (RMST.AND.CHRMMT) THEN
1565:              CALL CHRMS(JP,RMSD)1564:              CALL CHRMS(JP,RMSD)
1566:              IF (DEBUG) WRITE(MYUNIT,'(A,F15.5)')'RMSD = ',RMSD1565:              IF (DEBUG) WRITE(MYUNIT,'(A,F15.5)')'RMSD = ',RMSD
1567:              IF (RMSD.LE.RMSLIMIT) CALL SAVERMS(JP,POTEL,RMSD)1566:              IF (RMSD.LE.RMSLIMIT) CALL SAVERMS(JP,POTEL,RMSD)
1568: !                NRMS=NRMS+11567: !                NRMS=NRMS+1


r31387/mc_gbh.F90 2016-10-26 09:30:14.972637992 +0100 r31386/mc_gbh.F90 2016-10-26 09:30:16.940664362 +0100
 15: !   You should have received a copy of the GNU General Public License 15: !   You should have received a copy of the GNU General Public License
 16: !   along with this program; if not, write to the Free Software     16: !   along with this program; if not, write to the Free Software    
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19: SUBROUTINE MC_GBH(NSTEPS) 19: SUBROUTINE MC_GBH(NSTEPS)
 20:   ! 20:   !
 21:   USE PORFUNCS 21:   USE PORFUNCS
 22:   USE COMMONS, ONLY : MYNODE, MYUNIT, DEBUG, NATOMSALLOC, NQ, & 22:   USE COMMONS, ONLY : MYNODE, MYUNIT, DEBUG, NATOMSALLOC, NQ, &
 23:        TSTART, COORDS, LABELS, ECONV, BOXCENTROIDT,RANDMULTIPERMT, & 23:        TSTART, COORDS, LABELS, ECONV, BOXCENTROIDT,RANDMULTIPERMT, &
 24:        NPAR_GBH, TEMP, RMS, TARGET, HIT, QALCST, QALCS_NBRHD, & 24:        NPAR_GBH, TEMP, RMS, TARGET, HIT, QALCST, QALCS_NBRHD, &
 25:        QALCSV, GBH_RESTART, GBH_NREJMAX, GBH_NAVOID, & 25:        QALCSV, GBH_RESTART, GBH_NREJMAX, GBH_NAVOID
 26:        QALCS_SURFT, QALCS_SURF_MODE 
 27:   ! 26:   !
 28:   IMPLICIT NONE 27:   IMPLICIT NONE
 29:   ! 28:   !
 30: #ifdef MPI  29: #ifdef MPI 
 31:   INCLUDE 'mpif.h' 30:   INCLUDE 'mpif.h'
 32:   INTEGER MPIERR 31:   INTEGER MPIERR
 33: #endif 32: #endif
 34:   ! 33:   !
 35:   INTEGER, INTENT(IN) :: NSTEPS 34:   INTEGER, INTENT(IN) :: NSTEPS
 36:   ! 35:   !
 37:   LOGICAL :: RESET, DUPLICATE 36:   LOGICAL :: RESET
 38:   INTEGER :: I, ITRIAL, ITERNS, BRUN, QDONE, NQTOT, IPROC, IPROCLO, & 37:   INTEGER :: I, ITRIAL, ITERNS, BRUN, QDONE, NQTOT, IPROC, IPROCLO, &
 39:        NREJSTREAK, L0(NATOMSALLOC), LDUM(NATOMSALLOC) 38:        NREJSTREAK, L0(NATOMSALLOC)
 40:   DOUBLE PRECISION :: TIME, SCREENC(3*NATOMSALLOC), POTEL, R, & 39:   DOUBLE PRECISION :: TIME, SCREENC(3*NATOMSALLOC), POTEL, R, &
 41:        POTEL_LIST(NPAR_GBH+1), DPRAND, X0(3*NATOMSALLOC), E0, & 40:        POTEL_LIST(NPAR_GBH+1), DPRAND, X0(3*NATOMSALLOC), E0
 42:        XDUM(3*NATOMSALLOC), EDUM 
 43:   DOUBLE PRECISION, ALLOCATABLE :: AVOIDLIST(:) 41:   DOUBLE PRECISION, ALLOCATABLE :: AVOIDLIST(:)
 44:   ! 42:   !
 45:   COMMON /MYPOT/ POTEL 43:   COMMON /MYPOT/ POTEL
 46:   ! 44:   !
 47:   WRITE(MYUNIT, '(A)')  'mc_gbh> Calculating initial energy' 45:   WRITE(MYUNIT, '(A)')  'mc_gbh> Calculating initial energy'
 48:   !WRITE(MYUNIT, *)  'mc_gbh> NATOMSALLOC=', NATOMSALLOC 46:   !WRITE(MYUNIT, *)  'mc_gbh> NATOMSALLOC=', NATOMSALLOC
 49:   CALL FLUSH(MYUNIT) 47:   CALL FLUSH(MYUNIT)
 50:   CALL QUENCH(.FALSE.,1,ITERNS,TIME,BRUN,QDONE,SCREENC) 48:   CALL QUENCH(.FALSE.,1,ITERNS,TIME,BRUN,QDONE,SCREENC)
 51:   NQ(1) = 0 49:   NQ(1) = 0
 52:   WRITE(MYUNIT,& 50:   WRITE(MYUNIT,&
 63:   WRITE(MYUNIT, '(A,I6)') & 61:   WRITE(MYUNIT, '(A,I6)') &
 64:        'mc_gbh> Starting GBH loop of length ',NSTEPS 62:        'mc_gbh> Starting GBH loop of length ',NSTEPS
 65:   ! 63:   !
 66:   X0(1:3*NATOMSALLOC) = COORDS(1:3*NATOMSALLOC,1) 64:   X0(1:3*NATOMSALLOC) = COORDS(1:3*NATOMSALLOC,1)
 67:   L0(1:NATOMSALLOC) = LABELS(1:NATOMSALLOC,1) 65:   L0(1:NATOMSALLOC) = LABELS(1:NATOMSALLOC,1)
 68:   E0 = POTEL 66:   E0 = POTEL
 69:   ! 67:   !
 70:   NREJSTREAK=0 68:   NREJSTREAK=0
 71:   gbh_loop: DO ITRIAL=1,NSTEPS 69:   gbh_loop: DO ITRIAL=1,NSTEPS
 72:      ! 70:      !
 73:      WRITE(MYUNIT, '(A,I10)') & 71:      IF(QALCSV) WRITE(MYUNIT, '(A,I10)') &
 74:           'mc_gbh> Start of iteration ', ITRIAL      72:           'mc_gbh> Starting iteration ', ITRIAL     
 75:      ! 73:      !
 76:      RESET=.FALSE. 74:      RESET=.FALSE.
 77:      IF(GBH_RESTART) THEN 75:      IF(GBH_RESTART) THEN
 78:         IF(NREJSTREAK >= GBH_NREJMAX) RESET=.TRUE. 76:         IF(NREJSTREAK >= GBH_NREJMAX) RESET=.TRUE.
 79:         IF(GBH_NAVOID>0.AND.(RESET.OR.NREJSTREAK==1)) THEN 77:         IF(GBH_NAVOID>0.AND.(RESET.OR.NREJSTREAK==1)) THEN
 80:            CALL PROCESS_AVOID_LIST(E0,GBH_NAVOID,AVOIDLIST,RESET) 78:            CALL PROCESS_AVOID_LIST(E0,GBH_NAVOID,AVOIDLIST,RESET)
 81:         ENDIF 79:         ENDIF
 82:         IF(RESET) THEN 80:         IF(RESET) THEN
 83:            E0=1.0D+99 ! Guarantee move by inflating Markov energy 81:            E0=1.0D+99 ! Guarantee move by inflating Markov energy
 84:            NREJSTREAK=0            82:            NREJSTREAK=0           
 95:      ! Quench perturbed state 93:      ! Quench perturbed state
 96:      NQ(1) = NQ(1) + 1 94:      NQ(1) = NQ(1) + 1
 97:      CALL QUENCH(.FALSE.,1,ITERNS,TIME,BRUN,QDONE,SCREENC) 95:      CALL QUENCH(.FALSE.,1,ITERNS,TIME,BRUN,QDONE,SCREENC)
 98:      WRITE(MYUNIT,& 96:      WRITE(MYUNIT,&
 99:           '(A,I10,A,G20.10,A,I5,A,G12.5,A,G20.10,A,F10.1)') & 97:           '(A,I10,A,G20.10,A,I5,A,G12.5,A,G20.10,A,F10.1)') &
100:           'Qu ',NQ(1),' E= ',POTEL,' steps= ',ITERNS,' RMS= ',& 98:           'Qu ',NQ(1),' E= ',POTEL,' steps= ',ITERNS,' RMS= ',&
101:           RMS,' E0=',E0,' t= ',TIME-TSTART 99:           RMS,' E0=',E0,' t= ',TIME-TSTART
102:      CALL FLUSH(MYUNIT)100:      CALL FLUSH(MYUNIT)
103:      !101:      !
104: #ifdef MPI102: #ifdef MPI
105:      IF (QALCS_SURF_MODE == -1) THEN 
106:         ! 
107:         ! Back-up randomly perturbd state on each PROC 
108:         XDUM(1:3*NATOMSALLOC) = COORDS(1:3*NATOMSALLOC,1) 
109:         LDUM(1:NATOMSALLOC) = LABELS(1:NATOMSALLOC,1) 
110:         EDUM = POTEL 
111:         ! 
112:         proc_loop: DO IPROC=0,NPAR_GBH-1 
113:            ! 
114:            WRITE(MYUNIT,'(A,I3)') & 
115:                 'mc_gbh> Refining state from IPROC= ',IPROC+1 
116:            ! 
117:            ! Reinstate the random state on IPROC 
118:            IF(MYNODE == IPROC) THEN 
119:               COORDS(1:3*NATOMSALLOC,1) = XDUM(1:3*NATOMSALLOC) 
120:               LABELS(1:NATOMSALLOC,1) = LDUM(1:NATOMSALLOC) 
121:               POTEL = EDUM 
122:            ENDIF 
123:            ! 
124:            ! Broadcast state from IPROC 
125:            CALL MPI_BCAST(COORDS(:,1),3*NATOMSALLOC,& 
126:                 MPI_DOUBLE_PRECISION,IPROC,MPI_COMM_WORLD,MPIERR) 
127:            CALL MPI_BCAST(LABELS(:,1),NATOMSALLOC,MPI_INTEGER,& 
128:                 IPROC,MPI_COMM_WORLD,MPIERR) 
129:            CALL MPI_BCAST(POTEL,1,MPI_DOUBLE_PRECISION,& 
130:                 IPROC,MPI_COMM_WORLD,MPIERR) 
131:            ! 
132:            ! Accumulate POTEL_LIST with randomly perturbed  
133:            ! energies and check for duplicates 
134:            POTEL_LIST(IPROC+1) = POTEL 
135:            DUPLICATE=.FALSE. 
136:            IF(IPROC>0) THEN 
137:               DUPLICATE=.FALSE. 
138:               check_dupes: DO I=IPROC,1,-1 
139:                  IF(ABS(POTEL - POTEL_LIST(I)) < ECONV) THEN 
140:                     DUPLICATE=.TRUE. 
141:                     WRITE(MYUNIT,'(A)') & 
142:                          'mc_gbh> Duplicate skipped.' 
143:                     EXIT check_dupes 
144:                  ENDIF 
145:               ENDDO check_dupes 
146:            ENDIF 
147:            ! 
148:            IF(.NOT.DUPLICATE .AND. & 
149:                 ABS(POTEL - E0) > ECONV) THEN 
150:               ! Do surface refinement (parallelised!) 
151:               CALL QALCS_SURF(1, ITERNS, TIME, BRUN, QDONE, SCREENC) 
152:               ! Update back-up on IPROC 
153:               IF(MYNODE == IPROC) THEN 
154:                  XDUM(1:3*NATOMSALLOC) = COORDS(1:3*NATOMSALLOC,1) 
155:                  LDUM(1:NATOMSALLOC) = LABELS(1:NATOMSALLOC,1) 
156:                  EDUM = POTEL 
157:               ENDIF 
158:            ELSE ! Just inflate energy to prohibit selection 
159:               IF(MYNODE == IPROC) EDUM = 1.0D+99  
160:            ENDIF 
161:            ! 
162:         ENDDO proc_loop 
163:         ! 
164:         ! Reset to back-up on each PROC 
165:         COORDS(1:3*NATOMSALLOC,1) = XDUM(1:3*NATOMSALLOC) 
166:         LABELS(1:NATOMSALLOC,1) = LDUM(1:NATOMSALLOC) 
167:         POTEL = EDUM 
168:         ! 
169:      ENDIF 
170:      ! 
171:      ! Gather all parallel quench energies on master node.103:      ! Gather all parallel quench energies on master node.
172:      CALL MPI_GATHER(POTEL, 1, MPI_DOUBLE_PRECISION, &104:      CALL MPI_GATHER(POTEL, 1, MPI_DOUBLE_PRECISION, &
173:           POTEL_LIST(1:NPAR_GBH), 1, MPI_DOUBLE_PRECISION, 0, &105:           POTEL_LIST(1:NPAR_GBH), 1, MPI_DOUBLE_PRECISION, 0, &
174:           MPI_COMM_WORLD, MPIERR)106:           MPI_COMM_WORLD, MPIERR)
175: #else107: #else
176:      POTEL_LIST(1) = POTEL108:      POTEL_LIST(1) = POTEL
177: #endif109: #endif
178:      !110:      !
179:      IF(MYNODE==0) THEN111:      IF(MYNODE==0) THEN
180:         !112:         !
182:         IF (TARGET) IPROCLO=MINLOC(POTEL_LIST(1:NPAR_GBH), 1)-1114:         IF (TARGET) IPROCLO=MINLOC(POTEL_LIST(1:NPAR_GBH), 1)-1
183:         !115:         !
184:         ! Choose an IPROC.116:         ! Choose an IPROC.
185:         CALL CHOOSE_FROM_LIST(NPAR_GBH+1,POTEL_LIST,IPROC)117:         CALL CHOOSE_FROM_LIST(NPAR_GBH+1,POTEL_LIST,IPROC)
186:         !118:         !
187:         IF(IPROC>NPAR_GBH) THEN ! Reject all119:         IF(IPROC>NPAR_GBH) THEN ! Reject all
188:            COORDS(1:3*NATOMSALLOC,1) = X0(1:3*NATOMSALLOC)120:            COORDS(1:3*NATOMSALLOC,1) = X0(1:3*NATOMSALLOC)
189:            LABELS(1:NATOMSALLOC,1) = L0(1:NATOMSALLOC) 121:            LABELS(1:NATOMSALLOC,1) = L0(1:NATOMSALLOC) 
190:            POTEL = E0122:            POTEL = E0
191:            IPROC = 0123:            IPROC = 0
192:            WRITE(MYUNIT, '(A)') 'mc_gbh> Rejected all trial states!'124:            !WRITE(MYUNIT, '(A)') 'mc_gbh> Trial rejected.'     
193:         ELSE125:         ELSE
194:            WRITE(MYUNIT, '(A,I3)') &126:            !WRITE(MYUNIT, '(A)') 'mc_gbh> Trial accepted.'    
195:                 'mc_gbh> Accepted IPROC= ', IPROC     
196:            IPROC=IPROC-1127:            IPROC=IPROC-1
197:         ENDIF128:         ENDIF
198:         !129:         !
199:      ENDIF130:      ENDIF
200:      !131:      !
201: #ifdef MPI 132: #ifdef MPI 
202:      ! Broadcast IPROC from master133:      ! Broadcast IPROC from master
203:      CALL MPI_BCAST(IPROC,1,MPI_INTEGER,&134:      CALL MPI_BCAST(IPROC,1,MPI_INTEGER,&
204:           0,MPI_COMM_WORLD,MPIERR)135:           0,MPI_COMM_WORLD,MPIERR)
205:      ! Broadcast state from IPROC136:      ! Broadcast state from IPROC
228:      ENDIF159:      ENDIF
229:      !160:      !
230:      IF(HIT) THEN161:      IF(HIT) THEN
231:         WRITE(MYUNIT,'(A,I3,A,I6)') &162:         WRITE(MYUNIT,'(A,I3,A,I6)') &
232:              'mc_gbh> Target hit stochastically by node ',IPROCLO+1,&163:              'mc_gbh> Target hit stochastically by node ',IPROCLO+1,&
233:              ' on trial ',ITRIAL164:              ' on trial ',ITRIAL
234:         EXIT gbh_loop165:         EXIT gbh_loop
235:      ENDIF166:      ENDIF
236:      !167:      !
237:      ! --- Deterministic refinement --------------------------168:      ! --- Deterministic refinement --------------------------
238:      IF(.NOT.HIT .AND. ABS(POTEL-E0)>ECONV) THEN169:      IF(.NOT.HIT .AND. QALCST .AND. QALCS_NBRHD>0 &
239:         IF(QALCS_SURFT .AND. QALCS_SURF_MODE /= -1) THEN170:           .AND. ABS(POTEL-E0)>ECONV) THEN
240:            CALL QALCS_SURF(1, ITERNS, TIME, BRUN, QDONE, SCREENC)171:         ! perform a variable neighbourhood search (in parallel)
241:            WRITE(MYUNIT,'(A,G20.10,A,F11.1)') &172:         CALL QALCS_PARALLEL(ITERNS, TIME, BRUN, QDONE, SCREENC)
242:            'mc_gbh> Refined E= ',POTEL,' t= ', TIME-TSTART 
243:         ENDIF 
244:         IF( QALCST .AND. QALCS_NBRHD>0) THEN 
245:            ! perform a variable neighbourhood search (in parallel) 
246:            CALL QALCS_PARALLEL(ITERNS, TIME, BRUN, QDONE, SCREENC) 
247:            WRITE(MYUNIT,'(A,G20.10,A,F11.1)') & 
248:            'mc_gbh> Biminimum E= ',POTEL,' t= ', TIME-TSTART 
249:         ENDIF 
250:      ENDIF173:      ENDIF
251:      !174:      !
252:      IF(HIT) EXIT gbh_loop175:      IF(HIT) EXIT gbh_loop
253:      !176:      !
254:      X0(1:3*NATOMSALLOC) = COORDS(1:3*NATOMSALLOC,1)177:      X0(1:3*NATOMSALLOC) = COORDS(1:3*NATOMSALLOC,1)
255:      L0(1:NATOMSALLOC) = LABELS(1:NATOMSALLOC,1)178:      L0(1:NATOMSALLOC) = LABELS(1:NATOMSALLOC,1)
256:      E0 = POTEL179:      E0 = POTEL
257:      !180:      !
258:   ENDDO gbh_loop181:   ENDDO gbh_loop
259:   !182:   !


r31387/QALCSearch.f90 2016-10-26 09:30:13.972624593 +0100 r31386/QALCSearch.f90 2016-10-26 09:30:15.916650641 +0100
193:      IF(QALCS_SYMT.OR.QALCS_SURFT) THEN193:      IF(QALCS_SYMT.OR.QALCS_SURFT) THEN
194:         E0 = POTEL194:         E0 = POTEL
195:         NQ0 = NQTOT195:         NQ0 = NQTOT
196:         IF(QALCS_SYMT) THEN ! Symmetrisation scheme196:         IF(QALCS_SYMT) THEN ! Symmetrisation scheme
197:            CALL QALCS_SYM(NP,ITER,TIME,BRUN,QDONE,SCREENC)197:            CALL QALCS_SYM(NP,ITER,TIME,BRUN,QDONE,SCREENC)
198:         ENDIF198:         ENDIF
199:         IF(QALCS_SURFT) THEN ! DLS-like surface optimisation.199:         IF(QALCS_SURFT) THEN ! DLS-like surface optimisation.
200:            CALL QALCS_SURF(NP,ITER,TIME,BRUN,QDONE,SCREENC)200:            CALL QALCS_SURF(NP,ITER,TIME,BRUN,QDONE,SCREENC)
201:         ENDIF201:         ENDIF
202:         IF( POTEL < E0 - ECONV .AND. QALCST) DONE = .FALSE.202:         IF( POTEL < E0 - ECONV .AND. QALCST) DONE = .FALSE.
 203:         CALL MYCPU_TIME(TIME)
 204:         WRITE(MYUNIT,'(A,F20.10,A,I9,A,F11.1)') &
 205:              'QALCS> Refined E= ',POTEL, &
 206:              ' after ',NQTOT-NQ0,' quenches t= ',TIME-TSTART        
203:      ENDIF207:      ENDIF
204:      !208:      !
205:   ENDDO209:   ENDDO
206:   !210:   !
207:   CALL MYCPU_TIME(TIME)211:   CALL MYCPU_TIME(TIME)
208:   WRITE(MYUNIT,'(A,F20.10,A,F11.1)') &212:   WRITE(MYUNIT,'(A,F20.10,A,F11.1)') &
209:        'QALCS> Final E= ',POTEL,' t= ',TIME-TSTART213:        'QALCS> Final E= ',POTEL,' t= ',TIME-TSTART
210:   !214:   !
211:   RETURN215:   RETURN
212:   !216:   !


r31387/QALCS_surface.F90 2016-10-26 09:30:13.472617894 +0100 r31386/QALCS_surface.F90 2016-10-26 09:30:15.212641207 +0100
  1: !   GMIN: A program for finding global minima  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/QALCS_surface.F90' in revision 31386
  2: !   Copyright (C) 1999-2006 David J. Wales  
  3: !   This file is part of GMIN. 
  4: ! 
  5: !   GMIN is free software; you can redistribute it and/or modIFy 
  6: !   it under the terms of the GNU General Public License as published by 
  7: !   the Free Software Foundation; either version 2 of the License, or 
  8: !   (at your option) any later version. 
  9: ! 
 10: !   GMIN is distributed in the hope that it will be useful, 
 11: !   but WITHOUT ANY WARRANTY; without even the implied warranty of 
 12: !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 13: !   GNU General Public License for more details. 
 14: ! 
 15: !   You should have received a copy of the GNU General Public License 
 16: !   along with this program; IF not, write to the Free Software 
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-130 US 
 18: ! 
 19: !============================================================= 
 20: !   All routines in this file were implemented by 
 21: !   Dmitri Schebarchov (ds656). 
 22: !=============================================================  
 23: ! 
 24: SUBROUTINE QALCS_SURF(NP,ITER,TIME,BRUN,QDONE,SCREENC) 
 25:   ! 
 26:   ! Quench-Assisted Local Combinatorial Search for surface vacancies. 
 27:   ! Atoms are systematically swapped with "vacancies" (type 0). 
 28:   ! 
 29:   USE COMMONS, ONLY : NATOMS, VSITES, NNLISTS, COORDS, & 
 30:        NQ, ECONV, MYUNIT, TSTART, QALCSV, BOXCENTROIDT, & 
 31:        NPAR_GBH, MYNODE, TARGET, HIT 
 32:   ! 
 33:   IMPLICIT NONE 
 34:   ! 
 35: #ifdef MPI 
 36:   INCLUDE 'mpif.h' 
 37:   INTEGER MPIERR 
 38: #endif 
 39:   ! 
 40:   ! Parse passed variables 
 41:   INTEGER, INTENT(IN) :: NP 
 42:   INTEGER, INTENT(INOUT) :: ITER, BRUN, QDONE ! for QUENCH 
 43:   DOUBLE PRECISION, INTENT(INOUT) :: TIME, SCREENC(3*NATOMS) ! for QUENCH 
 44:   ! 
 45:   LOGICAL ::  COMPLETE 
 46:   INTEGER :: I1,I2,I3,I33,K,NQTOT,ATOMS_SORTED_BY_NN(0:NATOMS), & 
 47:        NN_LOWEST, NVSITES, VSITE_WEIGHT, NAVPAIRS, NPPN, & 
 48:        MYI1, MYI2, IPROC, NQ0 
 49:   DOUBLE PRECISION :: X(3*NATOMS),E,POTEL,NN_AVE, NNDIST(3), & 
 50:        XBEST(3*NATOMS), EBEST, POTEL_LIST(NPAR_GBH) 
 51:   ! 
 52:   INTEGER, ALLOCATABLE :: AVPAIRS(:,:) 
 53:   ! 
 54:   ! Energy of COORDS from last quench. Common block in QUENCH. 
 55:   COMMON /MYPOT/ POTEL 
 56:   ! Total quench count. Commom block in MC. 
 57:   COMMON /TOT/ NQTOT 
 58:   ! 
 59:   X(1:3*NATOMS) = COORDS(1:3*NATOMS,NP) 
 60:   XBEST(1:3*NATOMS) = COORDS(1:3*NATOMS,NP) 
 61:   E = POTEL 
 62:   EBEST = POTEL 
 63:   COMPLETE = .FALSE. 
 64:   NQ0 = NQ(NP) 
 65:   ! 
 66:   WRITE(MYUNIT,'(A,F20.10,A,F11.1)') & 
 67:        'QALCS_surf> Initial E= ',POTEL,' t= ',TIME-TSTART 
 68:   ! 
 69:   DO WHILE (.NOT. COMPLETE) 
 70:      ! 
 71:      COMPLETE = .TRUE. 
 72:      ! 
 73:      CALL BUILD_NNLISTS(NP, ATOMS_SORTED_BY_NN, NNDIST) 
 74:      CALL GEN_VSITES(NP, ATOMS_SORTED_BY_NN, NNDIST, & 
 75:           NVSITES, VSITE_WEIGHT)  
 76:      ! 
 77:      ! Consider only atoms of coordination 6 and below 
 78:      ! NN_AVE is apparently too large as an upper bound...  
 79:      ! Perhaps try NN_LOWEST? (By analogy with MAX_WEIGHT for voids.) 
 80:      NN_LOWEST = NNLISTS(ATOMS_SORTED_BY_NN(1),0) 
 81:      !write(*,*) NN_LOWEST 
 82:      IF(NN_LOWEST <= VSITE_WEIGHT) THEN 
 83:      !IF(NN_LOWEST <= VSITE_WEIGHT + 1) THEN ! maybe better?!? 
 84:         REDUCE: DO I2=1,ATOMS_SORTED_BY_NN(0)  
 85:            I3=ATOMS_SORTED_BY_NN(I2) ! actual atom index 
 86:            IF(NNLISTS(I3,0) > NN_LOWEST) THEN 
 87:            !IF(NNLISTS(I3,0) > MAX(NN_LOWEST,6)) THEN 
 88:               ATOMS_SORTED_BY_NN(0) = I2-1 
 89:               EXIT REDUCE 
 90:            ENDIF 
 91:         ENDDO REDUCE 
 92:      ELSE 
 93:         ATOMS_SORTED_BY_NN(0) = 0 
 94:      ENDIF 
 95:      ! 
 96:      NAVPAIRS=ATOMS_SORTED_BY_NN(0)*NVSITES 
 97:      WRITE(MYUNIT,'(A,I3,A,I3,A,I6)') & 
 98:           'QALCS_surf> ',ATOMS_SORTED_BY_NN(0), & 
 99:           ' atom(s) and ',NVSITES,' v-site(s) => NAVPAIRS= ', & 
100:           NAVPAIRS 
101:      ! 
102:      IF(ALLOCATED(AVPAIRS)) DEALLOCATE(AVPAIRS) 
103:      ALLOCATE(AVPAIRS(NAVPAIRS,2)) 
104:      K=0 
105:      DO I1=1,NVSITES ! loop thru vacancies 
106:         DO I2=1,ATOMS_SORTED_BY_NN(0) ! loop thru atoms            
107:            I3=ATOMS_SORTED_BY_NN(I2) ! actual atom index 
108:            K = K+1 
109:            AVPAIRS(K,1:2) = (/I3,I1/) 
110:         ENDDO 
111:      ENDDO 
112:      IF(K /= NAVPAIRS) THEN 
113:         WRITE(MYUNIT,'(A)') 'QALCS_surf> ERROR: K /= NAVPAIRS' 
114:         STOP 
115:      ENDIF 
116:      ! 
117:      IF(MOD(NAVPAIRS,NPAR_GBH)==0) THEN 
118:         NPPN = NAVPAIRS/NPAR_GBH 
119:      ELSE 
120:         NPPN = NAVPAIRS/NPAR_GBH + 1 ! integer division!!! 
121:      ENDIF 
122:      MYI1 = MYNODE*NPPN+1 
123:      MYI2 = MIN(MYI1-1+NPPN, NAVPAIRS) 
124:      ! 
125:      DO I2=MYI1,MYI2 
126:         ! 
127:         I1=AVPAIRS(I2,2) ! vacancy index        
128:         I3=AVPAIRS(I2,1) ! atom index 
129:         ! 
130:         I33=3*(I3-1) 
131:         DO K=1,3 
132:            COORDS(I33+K,NP) = VSITES(I1,K) 
133:         ENDDO 
134:         ! 
135:         IF(BOXCENTROIDT) CALL BOXCENTROID(COORDS(:,NP)) 
136:         ! 
137:         NQTOT = NQTOT + 1 
138:         NQ(NP) = NQ(NP) + 1 
139:         CALL QUENCH(.FALSE.,NP,ITER,TIME,BRUN,QDONE,SCREENC) 
140:         IF(QALCSV) CALL PRINT_QUENCH(NP, ITER, '  ') 
141:         ! 
142:         IF(POTEL < EBEST - ECONV) THEN ! improvemenet 
143:            EBEST = POTEL 
144:            XBEST(1:3*NATOMS) = COORDS(1:3*NATOMS,NP) 
145:            WRITE(MYUNIT,'(A,F20.10,A,2(1X,I3))') & 
146:                 'QALCS_surf> New E_BEST= ', EBEST,& 
147:                 ' for A,V= ',I3,I1  
148:         ENDIF 
149:         ! 
150:         POTEL = E 
151:         COORDS(1:3*NATOMS,NP) = X(1:3*NATOMS) 
152:         ! 
153:      END DO 
154:      ! 
155:      IF(ALLOCATED(AVPAIRS)) DEALLOCATE(AVPAIRS) 
156:      ! 
157: #ifdef MPI 
158:      ! Gather all parallel quench energies on master node. 
159:      CALL MPI_GATHER(EBEST, 1, MPI_DOUBLE_PRECISION, & 
160:           POTEL_LIST, 1, MPI_DOUBLE_PRECISION, 0, & 
161:           MPI_COMM_WORLD, MPIERR) 
162:      IF(MYNODE==0) IPROC=MINLOC(POTEL_LIST, 1)-1 
163:      ! 
164:      ! Broadcast IPROC from master 
165:      CALL MPI_BCAST(IPROC,1,MPI_INTEGER,& 
166:           0,MPI_COMM_WORLD,MPIERR) 
167:      ! Broadcast state from IPROC 
168:      CALL MPI_BCAST(XBEST,3*NATOMS,& 
169:           MPI_DOUBLE_PRECISION,IPROC,MPI_COMM_WORLD,MPIERR) 
170:      CALL MPI_BCAST(EBEST,1,MPI_DOUBLE_PRECISION,& 
171:           IPROC,MPI_COMM_WORLD,MPIERR) 
172:      IF(TARGET) THEN 
173:         CALL MPI_BCAST(HIT,1,MPI_LOGICAL,& 
174:              IPROC,MPI_COMM_WORLD,MPIERR) 
175:         IF(HIT) WRITE(MYUNIT,'(A,I5)') & 
176:              'QALCS_surface> Target hit deterministically by node', & 
177:              IPROC+1 
178:      ENDIF 
179: #endif 
180:      IF(EBEST < E - ECONV) COMPLETE = .FALSE. 
181:      ! 
182:      POTEL = EBEST 
183:      COORDS(1:3*NATOMS,NP) = XBEST(1:3*NATOMS) 
184:      E = EBEST 
185:      X(1:3*NATOMS) = XBEST(1:3*NATOMS) 
186:      ! 
187:   END DO 
188:   ! 
189:   CALL MYCPU_TIME(TIME) 
190:   WRITE(MYUNIT,'(A,F20.10,A,I9,A,F11.1)') & 
191:        'QALCS_surf> Refined E= ',POTEL, & 
192:        ' after ',NQ(NP)-NQ0,' quenches, t= ',TIME-TSTART         
193:   ! 
194:   !WRITE(MYUNIT,'(A,F20.10,A,F11.1)') & 
195:   !     'QALCS_surf> Final E= ',POTEL,' t= ',TIME-TSTART   
196:   ! 
197:   RETURN 
198:   ! 
199: END SUBROUTINE QALCS_SURF 
200: ! 
201: !================================================================= 
202: ! 
203: SUBROUTINE BUILD_NNLISTS(NP, LIST, NNDIST) 
204:   ! 
205:   ! This non-standard way of builiding NNLISTS is tailored 
206:   ! for generating vacancies on a surface. One noteworthy side- 
207:   ! effect is that, for a pair of atoms (I and J), it can happen  
208:   ! that I is listed as a neighbour of J, but J is not listed as a  
209:   ! neighbour of I. 
210:   ! 
211:   ! Although the procedure does not distinguish between different 
212:   ! atomic species/types, it does respect atom groups: 
213:   ! the nearest-neighboour analysis is restricted to atoms in groups 
214:   ! 1 and 2; atoms in group 3 are completely ignored. Further, 
215:   ! only group-1 atoms will be listed in LIST on output; and  
216:   ! the minimum, average and maximum NN distances stored in NNDIST 
217:   ! will be computed exclusively from the local environment of  
218:   ! atoms in group 1. 
219:   ! 
220:   USE COMMONS, ONLY : NNLISTS,ATOMLISTS,INVATOMLISTS,NATOMS, &  
221:        NSPECIES, COORDS, MYUNIT, NNLIST_MAX 
222:   ! 
223:   IMPLICIT NONE 
224:   ! 
225:   INTEGER, INTENT(IN) :: NP 
226:   INTEGER, INTENT(OUT) :: LIST(0:NATOMS) 
227:   DOUBLE PRECISION, INTENT(OUT) :: NNDIST(3) 
228:   ! 
229:   LOGICAL :: GROWN, GAP, SELF 
230:   INTEGER :: T1,T2,G1,G2,J1,J2,I1,I2,I13,I23,K,L,M,N, IT 
231:   INTEGER :: LJ1,UJ1,LJ2,UJ2,G2START,IPAIR1(2),IPAIR2(2),NN_SURF 
232:   DOUBLE PRECISION :: DUMMY,DIST2,GAPTOL,RNGTOL,NN_AVE, & 
233:        NNDISTS(NATOMS,NNLIST_MAX) 
234:    
235:   ! 
236:   NNLISTS(1:NATOMS,0:NNLIST_MAX) = 0 
237:   NNDISTS(1:NATOMS,1:NNLIST_MAX) = 0.0D0 
238:   ! 
239:   ! Set the tolerance for defining a gap in the distribution of 
240:   ! nearest-neighbour distances. GAPTOL = 1.1 
241:   ! maximum deviation from minimal NN distance: RNGTOL=1.457~((1+SQRT(2))/2)^2 
242:   !GAPTOL = 1.457D0 ! 15% jump? 
243:   !RNGTOL = 2.0D0 ! 1+0.75*fccgap 
244:   GAPTOL = 1.333D0 
245:   RNGTOL = 1.457d0 
246:   NNDIST(1) = 1.0D9 ! initialise minimal value squared NN distance 
247:   ! 
248:   ! 1st pass: fill and sort lists 
249:   ! 
250:   DO T1=1,NSPECIES(0) ! Double loop over ALL types 
251:      DO T2=1,T1 
252:         ! 
253:         DO G1=1,2 ! Double loop over groups 1 & 2 (not 3) 
254:            ! 
255:            IF(T1 == T2) THEN 
256:               G2START = G1 
257:            ELSE 
258:               G2START = 1 
259:            ENDIF 
260:            ! 
261:            DO G2 = G2START,2 
262:               ! 
263:               ! Set bounds for outer loop over atoms 
264:               LJ1=1 
265:               UJ1=ATOMLISTS(T1,G1,0) 
266:               IF(T1==T2.AND.G1==G2) THEN 
267:                  SELF = .TRUE. 
268:                  UJ1 = UJ1 - 1 ! Exclude distance to itself 
269:               ELSE 
270:                  SELF = .FALSE. 
271:               ENDIF 
272:               ! 
273:               DO J1=LJ1,UJ1 ! Double loop over atoms 
274:                  I1=ATOMLISTS(T1,G1,J1) 
275:                  I13 = 3*(I1-1) 
276:                  ! 
277:                  ! Set bounds for inner loop over atoms 
278:                  IF(SELF) THEN 
279:                     LJ2=J1+1 
280:                     UJ2=UJ1+1 
281:                  ELSE 
282:                     LJ2=1 
283:                     UJ2=ATOMLISTS(T2,G2,0) 
284:                  ENDIF 
285:                  ! 
286:                  DO J2=LJ2,UJ2 
287:                     I2=ATOMLISTS(T2,G2,J2) 
288:                     I23 = 3*(I2-1) 
289:                     ! 
290:                     IPAIR1=(/I1,I2/) ! pair to be swapped 
291:                     IPAIR2=(/I2,I1/) ! useful later 
292:                     ! 
293:                     DIST2=0.0D0 
294:                     DO K=1,3 
295:                        DIST2 = DIST2 + & 
296:                             (COORDS(I23+K,NP)-COORDS(I13+K,NP))**2 
297:                     ENDDO 
298:                     IF(DIST2<NNDIST(1)) NNDIST(1) = DIST2 
299:                     ! 
300:                     DO K=1,2 ! loop through the pair of lists 
301:                        ! Check if list is still growing 
302:                        IF(NNLISTS(IPAIR1(K),0) < NNLIST_MAX) THEN 
303:                           NNLISTS(IPAIR1(K),0) = & 
304:                                NNLISTS(IPAIR1(K),0) + 1 
305:                           GROWN = .TRUE. 
306:                        ELSEIF( DIST2 < NNDISTS( IPAIR1(K), & 
307:                             NNLISTS(IPAIR1(K),0) ) ) THEN 
308:                           GROWN = .TRUE. 
309:                        ELSE 
310:                           GROWN = .FALSE. 
311:                        ENDIF 
312:                        ! Update the last entry if required 
313:                        IF( GROWN ) THEN 
314:                           M = NNLISTS(IPAIR1(K),0) 
315:                           NNLISTS(IPAIR1(K),M) = IPAIR2(K) 
316:                           NNDISTS(IPAIR1(K),M) = DIST2 
317:                           ! 
318:                           !ds656> testing... 
319:                           !write(MYUNIT,*) & 
320:                           !     'build_nnlists> list for atom',& 
321:                           !     IPAIR1(K),' : ',& 
322:                           !     (NNLISTS(IPAIR1(K),IT), & 
323:                           !     IT=1,NNLISTS(IPAIR1(K),0) ) 
324:                           !write(MYUNIT,*) & 
325:                           !     'build_nnlists> dists: ',& 
326:                           !     (NNDISTS(IPAIR1(K),IT), & 
327:                           !     IT=1,NNLISTS(IPAIR1(K),0) ) 
328:                           !<ds656 ...testing 
329:                           ! Shuffle the last entry if possible 
330:                           SHUFFLE: DO L=NNLISTS(IPAIR1(K),0),2,-1 
331:                              M=L-1 
332:                              IF(NNDISTS(IPAIR1(K),L) < & 
333:                                   NNDISTS(IPAIR1(K),M)) THEN 
334:                                 ! Swap atom indices in list 
335:                                 N = NNLISTS(IPAIR1(K),M) 
336:                                 NNLISTS(IPAIR1(K),M) = & 
337:                                      NNLISTS(IPAIR1(K),L)  
338:                                 NNLISTS(IPAIR1(K),L) = N 
339:                                 ! Swap distances 
340:                                 DUMMY = NNDISTS(IPAIR1(K),M) 
341:                                 NNDISTS(IPAIR1(K),M) = & 
342:                                      NNDISTS(IPAIR1(K),L)  
343:                                 NNDISTS(IPAIR1(K),L) = DUMMY 
344:                              ELSE 
345:                                 EXIT SHUFFLE 
346:                              ENDIF 
347:                           ENDDO SHUFFLE 
348:                           ! 
349:                        ENDIF ! check for change in list 
350:                        ! 
351:                     ENDDO ! Looped over two lists 
352:                     ! 
353:                  ENDDO ! Double loop over atom indices 
354:               ENDDO 
355:               ! 
356:            ENDDO ! Double loop over groups 
357:         ENDDO 
358:         ! 
359:      ENDDO ! Double loop over types 
360:   ENDDO 
361:   ! 
362:   !ds656> testing... 
363:   !WRITE(MYUNIT,*) 'gen_nnlists> NN counts after 1st pass:', & 
364:   !     (NNLISTS(J1,0), J1 = 1, NATOMS) 
365:   !DO J1=1,NATOMS 
366:   !   write(MYUNIT,'(A,I3,A)',ADVANCE='NO') & 
367:   !        'build_nnlists> List for ', J1,':' 
368:   !   DO J2=1,NNLISTS(J1,0) 
369:   !      write(MYUNIT,'(I3,A,F6.4,A)',ADVANCE='NO')  & 
370:   !           NNLISTS(J1,J2), '(', NNDISTS(J1,J2),')' 
371:   !   ENDDO 
372:   !   WRITE(MYUNIT,*) 
373:   !ENDDO 
374:   !<ds656 ...testing 
375:   ! 
376:   ! 2nd pass: remove comparatively distant neighbours from lists  
377:   L=0 
378:   M=0 
379:   LIST(0:NATOMS) = 0 ! initialise the sorted list of atom labels 
380:   NNDIST(2:3) = 0.0D0 ! initialise average and maximium values 
381:   DO J1=1,NATOMS 
382:      IF(INVATOMLISTS(J1,2) == 1) THEN ! consider group 1 only 
383:         ! Check that the current atom's nearest neighbour  
384:         ! distance is not much greater than the overall  
385:         ! smallest neighbour distance.  
386:         IF(NNDISTS(J1,1)/NNDIST(1) < RNGTOL) THEN 
387:            GAP=.FALSE. 
388:            DUMMY = NNDISTS(J1,1) 
389:            GAPSEARCH: DO J2=2,NNLISTS(J1,0) 
390:               ! Identify a gap by spanning the ordered list of  
391:               ! NNDISTs (*squared* distances) and locating the  
392:               ! first gap greater than GAPTOL 
393:               ! [ Note: RNGTOL = 1.457 ~ ((1+SQRT(2))/2)^2 ] 
394:               !IF(NNDISTS(J1,J2)/NNDISTS(J1,J2-1) < GAPTOL .AND. & 
395:               !     NNDISTS(J1,J2)/NNDISTS(J1,1) < RNGTOL) THEN    
396:               IF(NNDISTS(J1,J2)/NNDISTS(J1,1) < RNGTOL) THEN      
397:                  DUMMY = DUMMY + NNDISTS(J1,J2) 
398:                  IF(NNDISTS(J1,J2) > NNDIST(3)) & 
399:                       NNDIST(3) = NNDISTS(J1,J2) 
400:               ELSE 
401:                  GAP=.TRUE. 
402:                  K = J2 
403:                  EXIT GAPSEARCH 
404:               ENDIF 
405:            ENDDO GAPSEARCH 
406:         ELSE 
407:            GAP=.TRUE. 
408:            K=1 
409:         ENDIF 
410:         IF(GAP) THEN 
411:            DO J2=K,NNLISTS(J1,0) 
412:               NNLISTS(J1,J2) = 0 
413:               NNDISTS(J1,J2) = 0.0D0 
414:            ENDDO 
415:            K = K - 1 
416:            ! Update count of neighbours for atom 
417:            NNLISTS(J1,0) = K 
418:         ELSE 
419:            K = NNLISTS(J1,0) 
420:         ENDIF 
421:         ! Accumulate average coordination 
422:         L = L + K 
423:         ! Accumulate average distance 
424:         NNDIST(2) = NNDIST(2) + DUMMY 
425:         ! Add atom index in LIST, sorted by coordination 
426:         ! from lowest to highest. 
427:         M = M + 1 
428:         LIST(M) = J1 
429:         ORDER: DO J2=M,2,-1 
430:            IF(NNLISTS(LIST(J2),0) < NNLISTS(LIST(J2-1),0)) THEN 
431:               K = LIST(J2) 
432:               LIST(J2) = LIST(J2-1) 
433:               LIST(J2-1) = K 
434:            ELSE 
435:               EXIT ORDER 
436:            ENDIF 
437:         ENDDO ORDER 
438:      !ELSE ! Set to zero and ignore in the future 
439:      !   NNLISTS(J1,0:NNLIST_MAX) = 0 
440:      !   NNDISTS(J1,1:NNLIST_MAX) = 0.0D0 
441:      ENDIF ! Group check 
442:   ENDDO 
443:   ! 
444:   NNDIST(2) = NNDIST(2)/DBLE(L) 
445:   NN_AVE = DBLE(L)/DBLE(M) 
446:   ! 
447:   !ds656> testing... 
448:   ! WRITE(MYUNIT,'(A)') 'build_nnlists> Final NN counts:' 
449:   ! DO K=1,NATOMS 
450:   !      WRITE(MYUNIT,'(I3)',ADVANCE='NO') NNLISTS(K,0) 
451:   ! ENDDO 
452:   ! WRITE(MYUNIT,'(A)') '' 
453:   !WRITE(MYUNIT,'(A)') 'gen_nnlists> Atoms sorted by NN count:' 
454:   !do J1=1,M ! M is the length of LIST set in the previous loop 
455:   !   j2= list(j1) 
456:   !   write(myunit,'(A,3(I5))') 'gen_nnlists> ',j1,j2, nnlists(j2,0) 
457:   !enddo 
458:   !<ds656 ...testing 
459:   ! 
460:   ! Choose cutoff for defining 'surface' atoms. 
461:   !NN_SURF = NINT(NN_AVE) 
462:   NN_SURF = 11 !?! Or NNLIST_MAX-1 ?!? 
463:   LIST(0) = M 
464:   SURFACE: DO K=1,M ! M <= NATOMS is the length of LIST 
465:      J1=LIST(K) ! Actual atom index 
466:      IF(NNLISTS(J1,0) > NN_SURF) THEN 
467:         LIST(0) = K-1 ! store count for 'surface' atoms 
468:         EXIT SURFACE  
469:      ENDIF       
470:   ENDDO SURFACE 
471:   ! 
472:   WRITE(MYUNIT,'(A,I3,1X,F6.3)') & 
473:        'build_nnlists> MIN and AVE NN-counts:', & 
474:        NNLISTS(LIST(1),0), NN_AVE 
475:   WRITE(MYUNIT,'(A,3(1X,F8.5))') & 
476:        'build_nnlists> MIN, AVE, MAX NN-distances: ', & 
477:        (SQRT(NNDIST(K)), K=1,3) 
478:   ! 
479:   RETURN 
480:   ! 
481: END SUBROUTINE BUILD_NNLISTS 
482: ! 
483: !================================================================= 
484: ! 
485: SUBROUTINE GEN_VSITES(NP, LIST, NNDIST, NVSITES, VSITE_WEIGHT) 
486:   ! 
487:   ! Generate vacant sites from NNLISTS 
488:   ! USE COMMONS, ONLY : NATOMS, NNLISTS, COORDS, NNLIST_MAX, & 
489:   !     VSITES, NVSITES_MAX, MYUNIT, ATOMLISTS, NSPECIES 
490:   USE COMMONS, ONLY : NATOMS, NNLISTS, COORDS, NNLIST_MAX, & 
491:       VSITES, NVSITES_MAX, MYUNIT, MIEFT, MIEF_NSITES, & 
492:       MIEF_SITES, MIEF_SIG, ATOMLISTS, NSPECIES 
493:   ! 
494:   IMPLICIT NONE 
495:   ! 
496:   INTEGER, INTENT(IN) :: NP 
497:   INTEGER, INTENT(IN) :: LIST(0:NATOMS) 
498:   DOUBLE PRECISION, INTENT(IN) :: NNDIST(3) 
499:   INTEGER, INTENT(OUT) :: NVSITES, VSITE_WEIGHT 
500:   INTEGER :: I0,I1,I2,I3,I13,I23,I33,J1,J2,J3,K, & 
501:        VSITE_WEIGHTS(NVSITES_MAX), NNLIST(0:NNLIST_MAX), & 
502:        VSITE_NNLIST(0:NNLIST_MAX), n0, n1, n2, n3 
503:   LOGICAL :: BAD, ADD 
504:   DOUBLE PRECISION :: VSITE(3), DUMMY, TOL1, TOL2, TOL3, TOL4 
505:   ! 
506:   VSITES(1:NVSITES_MAX,1:3) = 0.0D0 
507:   NVSITES=0 
508:   VSITE_WEIGHTS(1:NVSITES_MAX) = 0 
509:   ! 
510:   ! Thresholds for cluster-cluster atoms 
511:   !TOL1=0.375d0*NNDIST(1) ! squared rad. of circumsphere of tetrahedron (for overalap) 
512:   TOL1=NNDIST(1)*(0.5d0*(1.0D0+DSQRT(3.0D0/8.0D0)))**2  
513:   TOL2=NNDIST(3) ! Coordination determined by the largest NN dist. 
514:   !TOL2=1.8D0*NNDIST(1) ! This may be better?!? 
515:   ! 
516:   ! Thresholds for cluster-substrate atoms 
517:   IF(MIEFT) THEN 
518:      TOL3=0.64D0*MINVAL(MIEF_SIG)**2 ! adequate?!? 
519:      TOL4=1.44D0*MAXVAL(MIEF_SIG)**2 
520:   ENDIF 
521:   ! 
522:   DO I0=1,LIST(0) ! loop over atoms in list 
523:      ! 
524:      I1=LIST(I0) ! actual atom index 
525:      I13=3*(I1-1) 
526:      ! 
527:      DO J1=1,NNLISTS(I1,0) ! max = NNLIST_MAX 
528:         I2=NNLISTS(I1,J1) 
529:         I23=3*(I2-1) 
530:         ! 
531:         !WRITE(MYUNIT,'(A)') & 
532:         !     'gen_sites> New candidate v-site' 
533:         ! 
534:         ! Generate candidate vsite 
535:         DO K=1,3 ! XYZ coordinates for potential vsite 
536:            VSITE(K) = & 
537:                 2.0D0*COORDS(I13+K,NP) - COORDS(I23+K,NP) 
538:         ENDDO 
539:         VSITE_NNLIST(0) = 1 
540:         VSITE_NNLIST(1) = I1 
541:         VSITE_NNLIST(2:NNLIST_MAX) = 0 
542:         ! 
543:         BAD=.FALSE. 
544:         ! 
545:         ! Check distance from vsite to other neighbouring atoms 
546:         BADNESS1: DO J2=1,NNLISTS(I1,0) 
547:            IF(J2 /= J1) THEN 
548:               I2=NNLISTS(I1,J2) 
549:               I23=3*(I2-1) 
550:               DUMMY = 0.0D0 
551:               DO K=1,3 
552:                  DUMMY = DUMMY + & 
553:                       (COORDS(I23+K,NP)-VSITE(K))**2 
554:               ENDDO 
555:               IF(DUMMY < TOL1) THEN ! overlap 
556:                  BAD = .TRUE. 
557:                  EXIT BADNESS1 
558:               ELSEIF(DUMMY < TOL2) THEN ! add as neighbour 
559:                  IF(VSITE_NNLIST(0) < NNLIST_MAX) THEN  
560:                     VSITE_NNLIST(0) = VSITE_NNLIST(0) + 1 
561:                     VSITE_NNLIST(VSITE_NNLIST(0)) = I2 
562:                     !WRITE(MYUNIT,'(A,I3,A,f10.5)') & 
563:                     !     'gen_sites> Added atom ',I2,' as nbro.',dsqrt(dummy) 
564:                  ELSE 
565:                     !WRITE(MYUNIT,'(A)') & 
566:                     !     'gen_sites> WARNING(1): exceeded max nbro count' 
567:                     !WRITE(MYUNIT,'(A,I3,A,f10.5)') & 
568:                     !     'gen_sites> Not added atom ',I2,' as nbro.',dsqrt(dummy) 
569:                  ENDIF 
570:                  ! -------------------------------------------------- 
571:                  ! Also consider nbros of I2 that are not nbros of I1 
572:                  ! -------------------------------------------------- 
573:                  CALL REDUCE_NNLIST(NNLIST_MAX, & 
574:                       NNLISTS(I1,0:NNLIST_MAX), & 
575:                       NNLISTS(I2,0:NNLIST_MAX), & 
576:                       NNLIST(0:NNLIST_MAX) ) 
577:                  DO J3=1,NNLIST(0) 
578:                     I3=NNLIST(J3) 
579:                     I33=3*(I3-1) 
580:                     DUMMY = 0.0D0 
581:                     DO K=1,3 
582:                        DUMMY = DUMMY + & 
583:                             (COORDS(I33+K,NP)-VSITE(K))**2 
584:                     ENDDO 
585:                     IF(DUMMY < TOL1) THEN ! overlap! 
586:                        ! This block turns out to be important, especially when 
587:                        ! distinction between 1st and 2nd nearest neighbours is 
588:                        ! ambiguous! 
589:                        BAD=.TRUE. 
590:                        EXIT BADNESS1 
591:                     ELSEIF( DUMMY < TOL2 ) THEN 
592:                        ADD=.TRUE. 
593:                        LK: DO K=1,VSITE_NNLIST(0) 
594:                           IF(VSITE_NNLIST(K) == I3) THEN 
595:                              ADD = .FALSE. 
596:                              EXIT LK 
597:                           ENDIF 
598:                        ENDDO LK 
599:                        IF(ADD) THEN  
600:                           IF (VSITE_NNLIST(0) < NNLIST_MAX) THEN 
601:                              VSITE_NNLIST(0) = VSITE_NNLIST(0) + 1 
602:                              VSITE_NNLIST(VSITE_NNLIST(0)) = I3 
603:                              !WRITE(MYUNIT,'(A,I3,A,f10.5)') & 
604:                              !     'gen_sites> Added atom ',I3,' as nbro.',dsqrt(dummy) 
605:                           !ELSE 
606:                              !WRITE(MYUNIT,'(A)') & 
607:                              !     'gen_sites> WARNING(2): exceeded max nbro count' 
608:                              !WRITE(MYUNIT,'(A,I3,A,f10.5)') & 
609:                              !     'gen_sites> Not added atom ',I3,' as nbro.',dsqrt(dummy) 
610:  
611:                           ENDIF 
612:                        ENDIF 
613:                     ENDIF 
614:                  ENDDO 
615:                  ! -------------------------------------------------- 
616:               ENDIF 
617:               ! 
618:            ENDIF 
619:         ENDDO BADNESS1 
620:         ! 
621:         ! Check distance from VSITE to Mie site(s) 
622:         VSITE_WEIGHT = 0 
623:         IF(MIEFT .AND. .NOT. BAD) THEN 
624:            BADNESS2: DO J2=1,MIEF_NSITES 
625:               DUMMY = 0.0D0 
626:               DO K=1,3 
627:                  DUMMY = DUMMY + & 
628:                       (MIEF_SITES(J2,K)-VSITE(K))**2 
629:               ENDDO 
630:               IF(DUMMY < TOL3) THEN ! overlap 
631:                  BAD = .TRUE. 
632:                  EXIT BADNESS2 
633:               ELSEIF(DUMMY < TOL4) THEN ! add extra weight  
634:                  VSITE_WEIGHT = VSITE_WEIGHT + 1 
635:               ENDIF 
636:            ENDDO BADNESS2 
637:         ENDIF 
638:         ! 
639:         IF(.NOT. BAD) THEN 
640:            IF(NVSITES < NVSITES_MAX) THEN 
641:               NVSITES=NVSITES+1 
642:               VSITES(NVSITES,1:3) = VSITE(1:3) 
643:               !test 
644:               !WRITE(MYUNIT,'(A,3(1X,F10.5))') 'gen_sites> NEW V-site:',& 
645:               !     VSITE(1:3)  
646:               VSITE_WEIGHTS(NVSITES) = VSITE_NNLIST(0) + VSITE_WEIGHT 
647:            ELSE 
648:               WRITE(MYUNIT,'(A,I6)') & 
649:                    'gen_sites> NVSITES exceeded ', NVSITES_MAX 
650:               STOP 
651:            ENDIF 
652:            IF(VSITE_WEIGHTS(NVSITES).GE.12) THEN 
653:               WRITE(MYUNIT,'(A,I3)') 'gen_vsites> ***WARNING*** v-site coordination',& 
654:                    VSITE_WEIGHTS(NVSITES) 
655:               ! ds656> testing .. 
656:               ! WRITE(MYUNIT,'(A)', advance='no') 'gen_vsites> v-site nbros:' 
657: !               do n0=1,vsite_nnlist(0) 
658: !                  WRITE(MYUNIT,'(1x,i3)', advance='no') vsite_nnlist(n0) 
659: !               enddo 
660: !               write(myunit,*) 
661:  
662: !               WRITE(MYUNIT,'(A,I3)') 'gen_vsites> v-site generated from nbro of', I1 
663: !               WRITE(MYUNIT,'(A,I3)') 'gen_vsites> the nbro is atom',nnlists(i1,j1) 
664: !               WRITE(MYUNIT,'(A)', advance='no') 'gen_vsites> all nbros:' 
665: !               do n0=1,nnlists(i1,0) 
666: !                  WRITE(MYUNIT,'(1x,i3)', advance='no') nnlists(i1,n0) 
667: !               enddo 
668: !               write(myunit,*) 
669:  
670: !               write(*,'(I5)') natoms+1 
671: !               write(*,'(A)') 'Atoms and vacancies before merging' 
672: !               DO n0=1,NSPECIES(0) !ds656> Should not exceed 10            
673: !                  DO n1=1,ATOMLISTS(n0,1,0) 
674: !                     n2=ATOMLISTS(n0,1,n1) 
675: !                     WRITE(*,'(I2,3(1X,F20.10))') n0,& 
676: !                          (COORDS(n3,NP), n3=3*(n2-1)+1,3*n2) 
677: !                  ENDDO 
678: !               ENDDO 
679: !               WRITE(*,'(A,3(1X,F20.10))') ' 0',(VSITE(J2), J2=1,3) 
680: !              stop 
681:               ! <ds656 .. testing 
682:            ENDIF 
683:         ENDIF 
684:         ! 
685:      ENDDO ! loop over neighbour list 
686:      ! 
687:   ENDDO ! loop over atoms 
688:   ! 
689:   !ds656> testing... 
690:   !WRITE(MYUNIT,'(A)') 'gen_vsites> unsorted Vsites:'   
691:   !DO I1=1,NVSITES 
692:   !   WRITE(MYUNIT,'(A,I6,I3,3(1X,F10.5))') 'gen_sites>', & 
693:   !        I1, VSITE_WEIGHTS(I1), VSITES(I1,1:3) 
694:   !ENDDO 
695:   !<ds656 ...testing 
696:   ! 
697:   !WRITE(MYUNIT,'(A,I6)') 'gen_sites> V-site candidates: ', NVSITES 
698:   ! 
699:   ! Sort VSITES by weight from biggest to smallest (CHECK!) 
700:   CALL SORT3_V2(NVSITES,NVSITES_MAX,VSITE_WEIGHTS,VSITES) 
701:   ! weed out vacancies whose avergage coordination is below  
702:   ! 
703:   !ds656> testing... 
704:   !WRITE(MYUNIT,'(A)') 'gen_vsites> Vsites sorted by occupancy:'   
705:   !DO I1=1,NVSITES 
706:   !   WRITE(MYUNIT,'(A,I6,I3,3(1X,F10.5))') 'gen_sites>', & 
707:   !        I1, VSITE_WEIGHTS(I1), VSITES(I1,1:3) 
708:   !ENDDO 
709:   write(*,'(I5)') natoms+nvsites 
710:   write(*,'(A)') 'Atoms and all 1st generation vacancies' 
711:   DO I0=1,NSPECIES(0) !ds656> Should not exceed 10            
712:      DO I1=1,ATOMLISTS(I0,1,0) 
713:         J1=ATOMLISTS(I0,1,I1) 
714:         WRITE(*,'(A,I1,3(1X,F20.10))') 'A',I0,& 
715:              (COORDS(J2,NP), J2=3*(J1-1)+1,3*J1) 
716:      ENDDO 
717:   ENDDO 
718:   DO I1=1,NVSITES 
719:      WRITE(*,'(I2,3(1X,F20.10))') VSITE_WEIGHTS(I1),& 
720:           (VSITES(I1, J2), J2=1,3) 
721:   ENDDO 
722:   STOP 
723:   !<ds656 ...testing 
724:   ! 
725:   I2=NVSITES 
726:   IF(VSITE_WEIGHTS(1) > 1) THEN ! or 1? 
727:      ! Weed out all but the largest WEIGHTS (highest coordination) 
728:      !I3=MIN(VSITE_WEIGHTS(1),5) ! May be better?!? 
729:      I3=VSITE_WEIGHTS(1)  
730:      WEED: DO I1=I2,1,-1 
731:         IF(VSITE_WEIGHTS(I1) < I3) THEN 
732:            NVSITES = NVSITES-1 
733:            VSITE_WEIGHTS(I1) = 0 
734:            VSITES(I1,1:3) = 0.0D0 
735:         ELSE 
736:            EXIT WEED 
737:         ENDIF 
738:      ENDDO WEED 
739:      !VSITE_WEIGHT = VSITE_WEIGHTS(1) 
740:      VSITE_WEIGHT = I3 
741:      !WRITE(MYUNIT,'(A,I2)') & 
742:      !     'gen_sites> Highest v-site weight: ',VSITE_WEIGHT 
743:   ELSE 
744:      NVSITES = 0 
745:   ENDIF 
746:   ! 
747:   ! ds656> testing .. 
748:   !if(vsite_weights(1) .ge. 12) then 
749:      ! write(*,'(I5)') natoms+nvsites 
750:      ! write(*,'(A)') 'Atoms and vacancies before merging' 
751:      ! DO I0=1,NSPECIES(0) !ds656> Should not exceed 10            
752:      !    DO I1=1,ATOMLISTS(I0,1,0) 
753:      !       J1=ATOMLISTS(I0,1,I1) 
754:      !       WRITE(*,'(I2,3(1X,F20.10))') I0,& 
755:      !            (COORDS(J2,NP), J2=3*(J1-1)+1,3*J1) 
756:      !    ENDDO 
757:      ! ENDDO 
758:      ! DO I1=1,NVSITES 
759:      !    WRITE(*,'(A,3(1X,F20.10))') ' 0',(VSITES(I1, J2), J2=1,3) 
760:      ! ENDDO 
761:      !stop 
762:   !endif 
763:   ! <ds656 .. testing 
764:   ! 
765:   I1 = 1  
766:   DO WHILE(I1 < NVSITES) 
767:      I2 = I1+1 
768:      DO WHILE(I2 <= NVSITES) 
769:         DUMMY = 0.0D0 
770:         DO K=1,3 
771:            DUMMY = DUMMY + & 
772:                 (VSITES(I1,K)-VSITES(I2,K))**2 
773:         ENDDO 
774:         IF(DUMMY < TOL1) THEN ! Absorb I2 into I2 and decrement NVSITES 
775:            VSITES(I1,:) = 0.5D0*( VSITES(I1,:) + VSITES(I2,:) ) 
776:            VSITES(I2,:) = VSITES(NVSITES,:) 
777:            VSITE_WEIGHTS(I2) = VSITE_WEIGHTS(NVSITES) 
778:            VSITES(NVSITES,:) = 0.0D0 
779:            VSITE_WEIGHTS(NVSITES) = 0 
780:            NVSITES = NVSITES - 1 
781:         ELSE ! update I2 
782:            I2 = I2 + 1 
783:         ENDIF 
784:      ENDDO 
785:      I1 = I1 + 1 
786:   ENDDO 
787:   ! 
788:   IF(NVSITES > 0) THEN 
789:      WRITE(MYUNIT,'(A,I3,A,I3)') & 
790:           'gen_sites> Generated ', NVSITES,' v-site(s) of weight',& 
791:           VSITE_WEIGHT 
792:   ELSE 
793:      WRITE(MYUNIT,'(A,I6)') & 
794:           'gen_sites> Failed to generate adequate v-sites!'      
795:   ENDIF 
796:  
797:   !ds656> testing .. 
798:   !if(vsite_weight == 12) then 
799:   ! write(*,'(I5)') natoms+nvsites 
800:   ! write(*,'(A)') 'Final atoms and vacancies after merging' 
801:   ! DO I0=1,NSPECIES(0) !ds656> Should not exceed 10            
802:   !    DO I1=1,ATOMLISTS(I0,1,0) 
803:   !       J1=ATOMLISTS(I0,1,I1) 
804:   !       WRITE(*,'(I2,3(1X,F20.10))') I0,& 
805:   !            (COORDS(J2,NP), J2=3*(J1-1)+1,3*J1) 
806:   !    ENDDO 
807:   ! ENDDO 
808:   ! DO I1=1,NVSITES 
809:   !    WRITE(*,'(A,3(1X,F20.10))') '0',(VSITES(I1, J2), J2=1,3) 
810:   ! ENDDO 
811:   ! stop 
812:   !endif 
813:   !<ds656 .. testing 
814:   ! 
815:   RETURN 
816:   ! 
817: END SUBROUTINE GEN_VSITES 
818: ! 
819: SUBROUTINE REDUCE_NNLIST(NNLIST_MAX,LISTA,LISTB,LISTC) 
820:   ! 
821:   ! Put in LISTC all entries from LISTB that are not in LISTA. 
822:   ! Lists are not assumed to be sorted, so the procedure is slow. 
823:   ! 
824:   IMPLICIT NONE 
825:   ! 
826:   INTEGER, INTENT(IN) :: NNLIST_MAX 
827:   INTEGER, INTENT(IN) :: LISTA(0:NNLIST_MAX), LISTB(0:NNLIST_MAX) 
828:   INTEGER, INTENT(OUT) :: LISTC(0:NNLIST_MAX) 
829:   INTEGER :: IA, IB  
830:   LOGICAL :: SHARED 
831:   ! 
832:   LISTC(0:NNLIST_MAX) = 0 
833:   ! 
834:   LB: DO IB=1,LISTB(0) 
835:      SHARED = .FALSE. 
836:      LA: DO IA = 1, LISTA(0) 
837:         IF(LISTB(IB) == LISTA(IA)) THEN 
838:            SHARED = .TRUE. 
839:            EXIT LA 
840:         ENDIF 
841:      ENDDO LA 
842:      IF(.NOT. SHARED) THEN 
843:         LISTC(0) = LISTC(0) + 1 
844:         LISTC(LISTC(0)) = LISTB(IB) 
845:      ENDIF 
846:   ENDDO LB 
847:   ! 
848:   RETURN 
849:   ! 
850: END SUBROUTINE REDUCE_NNLIST 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0