hdiff output

r31975/gay-berne.f90 2017-02-27 09:30:11.736202307 +0000 r31974/gay-berne.f90 2017-02-27 09:30:11.960205241 +0000
4142:         PYA1bin(:,:)=0.5D04142:         PYA1bin(:,:)=0.5D0
4143:         DO J1=1,NATOMS/24143:         DO J1=1,NATOMS/2
4144:          RCUTARRAY(J1) = 2.0D04144:          RCUTARRAY(J1) = 2.0D0
4145:         END DO4145:         END DO
4146:          RCUT = 2.0D04146:          RCUT = 2.0D0
4147:       ENDIF4147:       ENDIF
4148: 4148: 
4149:       REALNATOMS = NATOMS/24149:       REALNATOMS = NATOMS/2
4150:       OFFSET     = 3 * REALNATOMS4150:       OFFSET     = 3 * REALNATOMS
4151: 4151: 
4152:         SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,NP)4152:                   SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,NP)
 4153:         y(:)=COORDS(:,NP)
4153: ! do the group moves if PERCGROUPT=.TRUE.4154: ! do the group moves if PERCGROUPT=.TRUE.
4154:    IF(PERCGROUPT) THEN4155:    IF(PERCGROUPT) THEN
4155:         nblocks=atomingroup(REALNATOMS) !the last number in the atomingroup array is the total number of groups as well4156:         nblocks=atomingroup(REALNATOMS) !the last number in the atomingroup array is the total number of groups as well
4156:         WRITE(MYUNIT,'(A,I6,A)') 'takestepelpsd> taking cooperative rotational and translational steps for ', nblocks, ' blocks' 4157:         WRITE(MYUNIT,*) 'takestepelpsd> taking cooperative rotational and translational steps for ', nblocks, ' blocks' 
4157:         if(allocated(rotationmatrix)) deallocate(rotationmatrix)4158:         if(allocated(rotationmatrix)) deallocate(rotationmatrix)
4158:         if(allocated(atomblockcentre)) deallocate(atomblockcentre)4159:         if(allocated(atomblockcentre)) deallocate(atomblockcentre)
4159:         allocate(atomblockcentre(nblocks,3),rotationmatrix(nblocks,3,3))4160:         allocate(atomblockcentre(nblocks,3),rotationmatrix(nblocks,3,3))
4160: !WRITE(MYUNIT,*) 'before rotation and translation:'4161: !WRITE(MYUNIT,*) 'before rotation and translation:'
4161: !WRITE(MYUNIT,*) y(:)4162: !WRITE(MYUNIT,*) y(:)
4162: OVRLPT=.TRUE. 
4163: 94  DO WHILE (OVRLPT) 
4164: 4163: 
4165:      y(:)=COORDS(:,NP)4164: 
4166: ! random rotation matrices for each group4165: ! random rotation matrices for each group
4167:      IF(ligrotscale>0.0D0) then4166:      IF(ligrotscale>0.0D0) then
4168:        DO i=1,nblocks  4167:        DO i=1,nblocks  
4169:         randomphi=(DPRAND()-0.5)*twopi*ligrotscale4168:         randomphi=(DPRAND()-0.5)*twopi*ligrotscale
4170:         randomtheta=(DPRAND()-0.5)*pi*ligrotscale4169:         randomtheta=(DPRAND()-0.5)*pi*ligrotscale
4171:         randompsi=(DPRAND()-0.5)*twopi*ligrotscale4170:         randompsi=(DPRAND()-0.5)*twopi*ligrotscale
4172:         st=sin(randomtheta)4171:         st=sin(randomtheta)
4173:         ct=cos(randomtheta)4172:         ct=cos(randomtheta)
4174:         sph=sin(randomphi)4173:         sph=sin(randomphi)
4175:         cph=cos(randomphi)4174:         cph=cos(randomphi)
4234:              randomz=DPRAND()4233:              randomz=DPRAND()
4235:              do k=1,REALNATOMS4234:              do k=1,REALNATOMS
4236:               if(i==atomingroup(k)) then4235:               if(i==atomingroup(k)) then
4237:                  y(3*k-2)=y(3*k-2)+2*randomx*ligtransstep*LOCALSTEP-ligtransstep*LOCALSTEP4236:                  y(3*k-2)=y(3*k-2)+2*randomx*ligtransstep*LOCALSTEP-ligtransstep*LOCALSTEP
4238:                  y(3*k-1)=y(3*k-1)+2*randomy*ligtransstep*LOCALSTEP-ligtransstep*LOCALSTEP4237:                  y(3*k-1)=y(3*k-1)+2*randomy*ligtransstep*LOCALSTEP-ligtransstep*LOCALSTEP
4239:                  y(3*k  )=y(3*k  )+2*randomz*ligtransstep*LOCALSTEP-ligtransstep*LOCALSTEP4238:                  y(3*k  )=y(3*k  )+2*randomz*ligtransstep*LOCALSTEP-ligtransstep*LOCALSTEP
4240:               end if4239:               end if
4241:              end do4240:              end do
4242:            end do4241:            end do
4243: 4242: 
4244:       COORDS(:,NP)=y(:) 
4245:  
4246: ! check for overlap 
4247:       DO J1 = 1, REALNATOMS 
4248:  
4249:          J3 = 3*J1 
4250:          J5 = OFFSET + J3 
4251:  
4252:             OVRLPT = .FALSE. 
4253:  
4254:             RI(:)  = COORDS(J3-2:J3,NP) 
4255:             P1(:)  = COORDS(J5-2:J5,NP) 
4256:  
4257: !     ROTATION MATRIX 
4258:  
4259:             CALL ELPSDRTMT (P1, PYA1bin(J1,:), AE) 
4260:  
4261:             DO J2 = 1, REALNATOMS 
4262:  
4263:                IF (J2 == J1) CYCLE 
4264:   
4265:                J4    = 3*J2 
4266:                J6    = OFFSET + J4 
4267:                RJ(:) = COORDS(J4-2:J4,NP) 
4268:                P2(:) = COORDS(J6-2:J6,NP) 
4269:  
4270: !     ROTATION MATRIX 
4271:  
4272:                CALL ELPSDRTMT (P2, PYA1bin(J2,:), BE) 
4273:  
4274:                RIJ    = RI - RJ 
4275:                ABSRIJ = DSQRT(DOT_PRODUCT(RIJ,RIJ)) 
4276:                          
4277:                IF (ABSRIJ < RCUT) THEN  
4278:  
4279: !     DETERMINE ELLIPTIC CONTACT FUNCTION 
4280:  
4281:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE, BE, RIJ, XMIN, FMIN) 
4282: 4243: 
4283:                   ECFVAL = - FMIN4244:       END IF ! PERCGROUPT
4284: ! allow for a slight overlap 4245:       COORDS(:,NP)=y(:)
4285:                   IF (ECFVAL < PYOVERLAPTHRESH) THEN 
4286:                      IF(DEBUG) WRITE(MYUNIT,*) 'takestepelpsd> PERCGROUP - particles overlapping', J1, J2, ECFVAL, ABSRIJ 
4287:                      OVRLPT = .TRUE. 
4288:                      GO TO 94 
4289:                   ENDIF 
4290:  
4291:                ENDIF 
4292:  
4293:             ENDDO  ! END LOOP OVER J2 
4294:       ENDDO  ! END LOOP OVER J1    
4295:  
4296:     END DO ! DO WHILE(OVERLAPT) 
4297:    END IF ! PERCGROUPT 
4298: !WRITE(MYUNIT,*) 'rotated and translated coordinates:'4246: !WRITE(MYUNIT,*) 'rotated and translated coordinates:'
4299: !WRITE(MYUNIT,*) y(:)4247: !WRITE(MYUNIT,*) y(:)
4300: 4248: 
4301:       DO J1 = 1,REALNATOMS4249:       DO J1 = 1,REALNATOMS
4302: 4250: 
4303:          J2     = 3*J14251:          J2     = 3*J1
4304:         IF (PARAMONOVPBCX) THEN4252:         IF (PARAMONOVPBCX) THEN
4305: !ensure x component of particle 1 vector is within BoxLx/2 of zero.4253: !ensure x component of particle 1 vector is within BoxLx/2 of zero.
4306: !if it isn't then subtract integer number of boxlx's such that it is.4254: !if it isn't then subtract integer number of boxlx's such that it is.
4307:                 COORDS(J2-2,NP)=COORDS(J2-2,NP)-BOXLX*NINT(COORDS(J2-2,NP)/BOXLX)4255:                 COORDS(J2-2,NP)=COORDS(J2-2,NP)-BOXLX*NINT(COORDS(J2-2,NP)/BOXLX)
4498:                         4446:                         
4499:                IF (ABSRIJ < RCUT) THEN 4447:                IF (ABSRIJ < RCUT) THEN 
4500: 4448: 
4501: !     DETERMINE ELLIPTIC CONTACT FUNCTION4449: !     DETERMINE ELLIPTIC CONTACT FUNCTION
4502: 4450: 
4503:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE, BE, RIJ, XMIN, FMIN)4451:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE, BE, RIJ, XMIN, FMIN)
4504: 4452: 
4505:                   ECFVAL = - FMIN4453:                   ECFVAL = - FMIN
4506: ! allow for a slight overlap 4454: ! allow for a slight overlap 
4507:                   IF (ECFVAL < PYOVERLAPTHRESH) THEN4455:                   IF (ECFVAL < PYOVERLAPTHRESH) THEN
4508:                         WRITE(MYUNIT,*) 'takestepelpsd> particles overlapping', J1, J2, ECFVAL, ABSRIJ4456:                         WRITE(MYUNIT,*) 'atoms overlapping', J1, J2, ECFVAL, ABSRIJ
4509:                      OVRLPT = .TRUE.4457:                      OVRLPT = .TRUE.
4510:                      GO TO 954458:                      GO TO 95
4511:                   ENDIF4459:                   ENDIF
4512: 4460: 
4513:                ENDIF4461:                ENDIF
4514: 4462: 
4515:             ENDDO  ! END LOOP OVER J24463:             ENDDO  ! END LOOP OVER J2
4516:          ENDDO  ! END WHILE4464:          ENDDO  ! END WHILE
4517:       ENDDO  ! END LOOP OVER J1   4465:       ENDDO  ! END LOOP OVER J1   
4518: 4466: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0