hdiff output

r31713/mc.F 2017-01-21 10:38:04.634750059 +0000 r31712/mc.F 2017-01-21 10:38:06.178885162 +0000
 22:       USE HINGE_MOVES 22:       USE HINGE_MOVES
 23:       use moves 23:       use moves
 24:       use grouprotmod 24:       use grouprotmod
 25:       use mbpolmod, only: mbpolstep 25:       use mbpolmod, only: mbpolstep
 26:  26: 
 27:       USE QMODULE , ONLY : QMIN, QMINP, INTEQMIN 27:       USE QMODULE , ONLY : QMIN, QMINP, INTEQMIN
 28:       USE modcharmm 28:       USE modcharmm
 29:       USE MODAMBER9, ONLY : MDSTEPT,CISARRAY1,CHIARRAY1,NOCISTRANSDNA,NOCISTRANSRNA, 29:       USE MODAMBER9, ONLY : MDSTEPT,CISARRAY1,CHIARRAY1,NOCISTRANSDNA,NOCISTRANSRNA,
 30:      &                      SETCHIRAL,AMCHPMAX,DOLIGMOVE,LIGMOVEFREQ, AMCHNMAX, ROTAMERT, 30:      &                      SETCHIRAL,AMCHPMAX,DOLIGMOVE,LIGMOVEFREQ, AMCHNMAX, ROTAMERT,
 31:      &                      E_IGB, E_BOND, E_ANGLE, E_DIHEDRAL, 31:      &                      E_IGB, E_BOND, E_ANGLE, E_DIHEDRAL,
 32:      &                      E_VDW, E_14_VDW, E_ELEC, E_14_ELEC, IGB, MACROIONT 32:      &                      E_VDW, E_14_VDW, E_ELEC, E_14_ELEC, IGB
 33:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, POT_ENE_REC_C 33:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, POT_ENE_REC_C
 34:       USE CHIRALITY, ONLY : INIT_CHIRAL, INIT_CIS_TRANS 34:       USE CHIRALITY, ONLY : INIT_CHIRAL, INIT_CIS_TRANS
 35:       USE porfuncs 35:       USE porfuncs
 36:       USE AMHGLOBALS, ONLY: NMRES,OMOVI,AVEP,NUMPRO,IRES 36:       USE AMHGLOBALS, ONLY: NMRES,OMOVI,AVEP,NUMPRO,IRES
 37:       USE AMH_INTERFACES, ONLY:E_WRITE 37:       USE AMH_INTERFACES, ONLY:E_WRITE
 38:       USE ROTAMER 38:       USE ROTAMER
 39:  39: 
 40:       IMPLICIT NONE 40:       IMPLICIT NONE
 41: #ifdef MPI 41: #ifdef MPI
 42:       INCLUDE 'mpif.h' 42:       INCLUDE 'mpif.h'
827:                ENDIF827:                ENDIF
828: 828: 
829: !829: !
830: ! csw34> Coordinates are saved so that moves can be undone830: ! csw34> Coordinates are saved so that moves can be undone
831: !831: !
832:                SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)832:                SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)
833: ! csw34> If you want to look at the effect of moves, you can dump out833: ! csw34> If you want to look at the effect of moves, you can dump out
834: ! the structure BEFORE the move here.834: ! the structure BEFORE the move here.
835: !                 CALL A9DUMPPDB(COORDS(:,JP),"beforemove")835: !                 CALL A9DUMPPDB(COORDS(:,JP),"beforemove")
836: !                 CALL CHARMMDUMP(COORDS(:,JP),'beforemove')836: !                 CALL CHARMMDUMP(COORDS(:,JP),'beforemove')
837:                IF (MACROIONT.AND..NOT.SANDBOXT) THEN837: 
838:                       CALL MACROION_MOVES(J1,COORDS(:,JP),movableatomlist,nmovableatoms,ligmovet,blockmovet,nblocks,atomsinblock,STEP(JP)) 
839:                ELSE 
840: !838: !
841: ! UNIVERSAL STEP TAKING839: ! UNIVERSAL STEP TAKING
842: !840: !
843: 841: 
844: ! csw34> Rigid body EXPAND moves842: ! csw34> Rigid body EXPAND moves
845:                      IF (RIGIDINIT.AND.EXPANDRIGIDT.AND.DOEXPANDRIGID) THEN843:                      IF (RIGIDINIT.AND.EXPANDRIGIDT.AND.DOEXPANDRIGID) THEN
846:                         CALL GENRIGID_EXPAND(COORDS(:,JP), EXPANDFACTOR) 844:                         CALL GENRIGID_EXPAND(COORDS(:,JP), EXPANDFACTOR) 
847:                      ENDIF845:                      ENDIF
848: ! csw34> Rigid body rotation moves846: ! csw34> Rigid body rotation moves
849:                      IF (RIGIDINIT.AND.ROTATERIGIDT.AND.DOROTATERIGID) THEN847:                      IF (RIGIDINIT.AND.ROTATERIGIDT.AND.DOROTATERIGID) THEN
850:                         CALL GENRIGID_ROTATE(COORDS(:,JP), ROTATEFACTOR)848:                         CALL GENRIGID_ROTATE(COORDS(:,JP), ROTATEFACTOR)
851:                      ENDIF849:                      ENDIF
852: ! mo361> Rigid body translation moves850: ! mo361> Rigid body translation moves
853:                      IF (RIGIDINIT.AND.TRANSLATERIGIDT.AND.DOTRANSLATERIGID) THEN851:                      IF (RIGIDINIT.AND.TRANSLATERIGIDT.AND.DOTRANSLATERIGID) THEN
854:                         CALL GENRIGID_TRANSLATE(COORDS(:,JP),TRANSLATEFACTOR)852:                         CALL GENRIGID_TRANSLATE(COORDS(:,JP),TRANSLATEFACTOR)
855:                      ENDIF853:                      ENDIF
856:                      IF (RIGIDINIT.AND.ROTATEHINGET.AND.DOROTATEHINGE) THEN854:                      IF (RIGIDINIT.AND.ROTATEHINGET.AND.DOROTATEHINGE) THEN
857:                         CALL HINGE_ROTATE(COORDS(:,JP), HINGEROTATEFACTOR)855:                         CALL HINGE_ROTATE(COORDS(:,JP), HINGEROTATEFACTOR)
858:                      ENDIF856:                      ENDIF
859:                END IF 
860: !857: !
861: ! CHARMM STEP TAKING858: ! CHARMM STEP TAKING
862: !859: !
863:                IF (CHRMMT) THEN860:                IF (CHRMMT) THEN
864:                   IF (CHMDT.AND.MOD(J1,CHMDFREQ).EQ.0) THEN861:                   IF (CHMDT.AND.MOD(J1,CHMDFREQ).EQ.0) THEN
865:                      CALL CHMD(JP)862:                      CALL CHMD(JP)
866:                   ELSE863:                   ELSE
867:                      IF (CHRIGIDTRANST.AND.MOD(J1,FTRANS).EQ.0) CALL MKRIGIDTRANS(JP)864:                      IF (CHRIGIDTRANST.AND.MOD(J1,FTRANS).EQ.0) CALL MKRIGIDTRANS(JP)
868:                      IF (CHRIGIDROTT.AND.MOD(J1,FROT).EQ.0) CALL MKRIGIDROT(JP)865:                      IF (CHRIGIDROTT.AND.MOD(J1,FROT).EQ.0) CALL MKRIGIDROT(JP)
869:                      IF (MOD(J1,CHFREQ).EQ.0) THEN866:                      IF (MOD(J1,CHFREQ).EQ.0) THEN
1125: !     csw34> Check distances for groups defined in movableatoms file1122: !     csw34> Check distances for groups defined in movableatoms file
1126: !            If A->B > ABTHRESH or A->C > ACTHRESH, the step is discarded and we try again 1123: !            If A->B > ABTHRESH or A->C > ACTHRESH, the step is discarded and we try again 
1127:                      IF (LOCALSAMPLET) THEN1124:                      IF (LOCALSAMPLET) THEN
1128:                         CALL A9DISTCHECK(COORDS(:,JP),DISTOK)1125:                         CALL A9DISTCHECK(COORDS(:,JP),DISTOK)
1129:                         IF (.NOT.DISTOK) COORDS(:,JP)=SAVECOORDS(:)1126:                         IF (.NOT.DISTOK) COORDS(:,JP)=SAVECOORDS(:)
1130:                      ELSE1127:                      ELSE
1131:                         DISTOK=.TRUE.1128:                         DISTOK=.TRUE.
1132:                      ENDIF1129:                      ENDIF
1133:                   ENDDO1130:                   ENDDO
1134: ! AMBER12 step taking1131: ! AMBER12 step taking
1135: ! with macroion model1132:                ELSE IF (AMBER12T) THEN
1136:                ELSE IF (AMBER12T.AND.MACROIONT) THEN 
1137:                   IF (MACROIONT) THEN 
1138:                       CALL MACROION_MOVES(J1,COORDS(:,JP),movableatomlist,nmovableatoms,ligmovet,blockmovet,nblocks,atomsinblock,STEP(JP)) 
1139:                   END IF 
1140: ! normal AMBER12 step taking (without MACROION model)               
1141:                ELSE IF (AMBER12T.AND..NOT.MACROIONT) THEN 
1142: ! Do group rotation moves1133: ! Do group rotation moves
1143:                   CALL CARTESIAN_SPHERE(COORDS(:,JP), STEP(JP))1134:                   CALL CARTESIAN_SPHERE(COORDS(:,JP), STEP(JP))
1144:                   IF (GROUPROTT.AND.DOGROUPROT) CALL GROUPROTSTEP(JP)1135:                   IF (GROUPROTT.AND.DOGROUPROT) CALL GROUPROTSTEP(JP)
1145:                   IF (DIHEDRALROTT.AND.DODIHEDRALROT) CALL DIHEDRALROTSTEP(DUMMY1,JP,DUMMY2,DUMMY3)1136:                   IF (DIHEDRALROTT.AND.DODIHEDRALROT) CALL DIHEDRALROTSTEP(DUMMY1,JP,DUMMY2,DUMMY3)
1146:                   IF (ROTAMER_MOVET) CALL ROTAMER_STEP(COORDS(:, JP))1137:                   IF (ROTAMER_MOVET) CALL ROTAMER_STEP(COORDS(:, JP))
1147: ! jdf43> MWFILM1138: ! jdf43> MWFILM
1148:                ELSE IF ((MWFILMT).AND.(PERIODIC)) THEN1139:                ELSE IF ((MWFILMT).AND.(PERIODIC)) THEN
1149:                   CALL MWSTEP(COORDS,JP,NPAR,NATOMS,STEP(JP),BOXLX,BOXLY,BOXLZ,LAT)1140:                   CALL MWSTEP(COORDS,JP,NPAR,NATOMS,STEP(JP),BOXLX,BOXLY,BOXLZ,LAT)
1150: ! vr274> DMACRYS steps taking1141: ! vr274> DMACRYS steps taking
1151:                ELSE IF (DMACRYST) THEN1142:                ELSE IF (DMACRYST) THEN


r31713/moves.f90 2017-01-21 10:38:04.910774289 +0000 r31712/moves.f90 2017-01-21 10:38:06.494912697 +0000
253:       XYZ(3 * I - 1) = (B * (W**2 + U**2) - V * (C*W + A*U - U*X - V*Y - W*Z)) * (1 - COS_THETA) + &253:       XYZ(3 * I - 1) = (B * (W**2 + U**2) - V * (C*W + A*U - U*X - V*Y - W*Z)) * (1 - COS_THETA) + &
254:                        Y * COS_THETA + &254:                        Y * COS_THETA + &
255:                        (-A*W + C*U + W*X - U*Z) * SIN_THETA255:                        (-A*W + C*U + W*X - U*Z) * SIN_THETA
256:       XYZ(3 * I    ) = (C * (U**2 + V**2) - W * (A*U + B*V - U*X - V*Y - W*Z)) * (1 - COS_THETA) + &256:       XYZ(3 * I    ) = (C * (U**2 + V**2) - W * (A*U + B*V - U*X - V*Y - W*Z)) * (1 - COS_THETA) + &
257:                        Z * COS_THETA + &257:                        Z * COS_THETA + &
258:                        (-B*U + A*V + U*Y - V*X) * SIN_THETA258:                        (-B*U + A*V + U*Y - V*X) * SIN_THETA
259:    END DO259:    END DO
260: 260: 
261: END SUBROUTINE ROTATION_ABOUT_AXIS261: END SUBROUTINE ROTATION_ABOUT_AXIS
262: 262: 
263: SUBROUTINE MACROION_MOVES(j1,y,movableatomlist1,nmovableatoms,ligmovet,blockmovet,nblocks,atomsinblock1,LOCALSTEP) 
264: ! sf344> this is a cooperative move set for the MACROION model. It could be adapted for other systems,  
265: ! where dynamically determined groups of atoms have to be moved together. 
266:  
267:    use modamber9 ! macroion related keywords are stored in this module 
268:    use commons, only : natoms, change_temp, newres_temp, perct, perccut, frozen 
269: !   use nblist, only : nbflag, skinnb 
270:    use porfuncs 
271:  
272:    implicit none 
273:  
274:  
275:    integer itime1, now(3), nmovableatoms,i,j,movableatomlist1(nmovableatoms),nblocks,atomsinblock1(nblocks),offset1,k,j1 
276:    double precision        :: y(3*natoms), grad(3*natoms), ereal,ligandcentre(3),ligandcoords(3),ligandcoordsrotated(3) 
277:    double precision        :: twopi,pi,DPRAND,randomphi,randompsi,randomtheta,DIST,dummyz,mindistance 
278:    double precision        :: st,ct,sph,cph,sps,cps, VECBAX, VECBAY, VECBAZ, randomx, randomy, randomz, dummyx, dummyy 
279:    double precision        :: distancematrix(natoms,natoms), ysave(3*natoms),cartstepsave,transstepsave,localstep 
280:    logical                 :: ligmovet,blockmovet,overlapt 
281:    logical                 :: includedatom(natoms),enabledmove(natoms) 
282:    character(len=10)       :: datechar,timechar,zonechar 
283:    integer                 :: values(8),iostatus,restemp,i1,currentresidue,resnumber,overlaparray(natoms,natoms) 
284:    integer                 :: overlaparraysave(natoms,natoms),nretries,loopcounter 
285:    character(len=10)       :: rotmaxchangestr,rotpselectstr,rotcutoffstr,rotcentrestr,rotoccuwstr 
286:    integer, allocatable    :: movableatomlist(:), atomsinblock(:)  
287:    double precision, allocatable :: atomblockcentre(:,:), rotationmatrix(:,:,:) 
288:  
289: ! sf344> for compatibility with dynamically determining macroion nearest neighbours and moving them with LIGMOVE 
290:    if(allocated(movableatomlist)) deallocate(movableatomlist) 
291:    allocate(movableatomlist(nmovableatoms)) 
292:    if(allocated(atomsinblock)) deallocate(atomsinblock) 
293:    allocate(atomsinblock(nblocks)) 
294:    if(allocated(rotationmatrix)) deallocate(rotationmatrix) 
295:    allocate(rotationmatrix(nblocks,3,3)) 
296:    if(allocated(atomblockcentre)) deallocate(atomblockcentre) 
297:    allocate(atomblockcentre(nblocks,3)) 
298:    movableatomlist(:)=movableatomlist1(:) 
299:    atomsinblock(:)=atomsinblock1(:) 
300:  
301:    cartstepsave=ligcartstep 
302:    transstepsave=ligtransstep 
303:  
304: ! parameters 
305:    pi=ATAN(1.0D0)*4 
306:    twopi=2.0D0*pi 
307:  
308: ! 1) RIGID BODY ROTATION 
309:  
310: ! sf344> work out for binary systems (macroion model) the bound counterions to each macroion, and move them cooperatively 
311: !  with BLOCKMOVE. 
312:  
313:  
314: ! Macroions are the first in the atom list. 
315:           nblocks=nmacroions 
316:           if(allocated(movableatomlist)) deallocate(movableatomlist) 
317:           if(allocated(atomsinblock)) deallocate(atomsinblock) 
318:           if(allocated(atomblockcentre)) deallocate(atomblockcentre) 
319:           if(allocated(rotationmatrix)) deallocate(rotationmatrix) 
320:           allocate(movableatomlist(natoms),atomsinblock(nblocks),atomblockcentre(nblocks,3),rotationmatrix(nblocks,3,3)) 
321:           movableatomlist(:)=0 
322:           atomsinblock(:)=0 
323:           atomblockcentre(:,:)=0.0D0 
324:           rotationmatrix(:,:,:)=0.0D0 
325:           offset1=1 
326:           includedatom(1:natoms)=.false. 
327:           do i=1,nmacroions 
328:                movableatomlist(offset1)=i 
329:                offset1=offset1+1 
330:                atomsinblock(i)=1 
331:              do j=nmacroions+1,natoms 
332:                distancematrix(i,j)=sqrt((y(3*i-2)-y(3*j-2))**2+(y(3*i-1)-y(3*j-1))**2+(y(3*i)-y(3*j))**2) 
333:                if(distancematrix(i,j)<macroiondist.and..not.includedatom(j))then 
334:                 includedatom(j)=.true. ! include one counterion in only one block! 
335:                 movableatomlist(offset1)=j 
336:                 offset1=offset1+1 
337:                 atomsinblock(i)=atomsinblock(i)+1 
338:                end if 
339:              end do 
340:           end do 
341: !          write(MYUNITNEW,*) 'dynamically determining bound counterions to each macroion and moving them together' 
342: !          write(MYUNITNEW,*) 'natoms, nblocks', natoms, nblocks 
343: !          write(MYUNITNEW,*) 'atomsinblock(:)', atomsinblock(:) 
344: !          write(MYUNITNEW,*) 'movableatomlist(:)', movableatomlist(:) 
345:         ysave(:)=y(:) 
346:  
347:    IF (LIGMOVET.AND.MOD(J1,ligmovefreq).EQ.0) THEN 
348:         overlapt=.true. 
349:         overlaparraysave(:,:)=0 
350:         overlaparray(:,:)=0 
351:         distancematrix(:,:)=0.0D0 
352:         ysave(:)=y(:) 
353:         IF(BLOCKMOVET) THEN 
354:          ! Try to avoid moves that cause extremely high initial forces (or cold fusion), 
355:          ! by redoing the rotation moves whenever two atoms get closer than 0.5 Angstroms. 
356:           do i=1,natoms 
357:            do j=1,natoms 
358:               if(i==j) cycle 
359:               distancematrix(i,j)=sqrt((y(3*i-2)-y(3*j-2))**2+(y(3*i-1)-y(3*j-1))**2+(y(3*i)-y(3*j))**2) 
360:               if((MACROIONT.and.distancematrix(i,j)<1.0D0).or.distancematrix(i,j)<0.7D0) then  
361:                 overlaparraysave(i,j)=1  !there's an initial overlap 
362: !                write(*,*) 'atoms close together before rotating ', i,j,distancematrix(i,j) 
363:               end if 
364:            end do 
365:           end do 
366:         END IF 
367:  
368:      nretries=0 
369:      loopcounter=0 
370:      overlapt=.true. 
371:      DO WHILE(OVERLAPT) 
372:        IF(nretries>10) then ! scale ligcartstep and ligtransstep 
373:           nretries=0 
374: !          ligcartstep=ligcartstep*0.5 
375: !          ligtransstep=ligtransstep*0.5 
376: !         WRITE(MYUNITNEW,*) 'scaling down ligcartstep and ligtransstep', ligcartstep, ligtransstep 
377:        END IF 
378:        IF(loopcounter>1000) exit !give up taking non-overlapping moves 
379:        overlapt=.false. 
380:        DO i=1,nblocks ! nblocks=1 if BLOCKMOVET=.FALSE.  
381:           randomphi=(DPRAND()-0.5)*twopi*ligrotscale*LOCALSTEP 
382:           randomtheta=(DPRAND()-0.5)*pi*ligrotscale*LOCALSTEP 
383:           randompsi=(DPRAND()-0.5)*twopi*ligrotscale*LOCALSTEP 
384:           st=sin(randomtheta) 
385:           ct=cos(randomtheta) 
386:           sph=sin(randomphi) 
387:           cph=cos(randomphi) 
388:           sps=sin(randompsi) 
389:           cps=cos(randompsi) 
390:  
391:           rotationmatrix(i,1,1)=cps*cph-ct*sph*sps 
392:           rotationmatrix(i,2,1)=cps*sph+ct*cph*sps 
393:           rotationmatrix(i,3,1)=sps*st 
394:           rotationmatrix(i,1,2)=-sps*cph-ct*sph*cps 
395:           rotationmatrix(i,2,2)=-sps*sph+ct*cph*cps 
396:           rotationmatrix(i,3,2)=cps*st 
397:           rotationmatrix(i,1,3)=st*sph 
398:           rotationmatrix(i,2,3)=-st*cph 
399:           rotationmatrix(i,3,3)=ct 
400: !        WRITE(*,*) 'rotation matrix: ', rotationmatrix(i,:,:) 
401:        END DO 
402:  
403: ! work out the centre of coordinates for the ligand 
404:        atomblockcentre(:,:)=0.0D0 
405:        offset1=0 
406: !      overlapt=.false. 
407:        do i=1,nblocks 
408: !          if(.not.enabledmove(i)) cycle 
409:           if(i>1) offset1=offset1+atomsinblock(i-1) 
410:           do k=1,atomsinblock(i) 
411:              j=movableatomlist(k+offset1) 
412: !            WRITE(*,*) k, j 
413:              atomblockcentre(i,1)=atomblockcentre(i,1)+y(3*j-2) 
414:              atomblockcentre(i,2)=atomblockcentre(i,2)+y(3*j-1) 
415:              atomblockcentre(i,3)=atomblockcentre(i,3)+y(3*j  ) 
416:           end do 
417:           atomblockcentre(i,:)=atomblockcentre(i,:)/atomsinblock(i) 
418:           do k=1,atomsinblock(i) 
419:              j=movableatomlist(k+offset1) 
420: ! move each block to the origin 
421:              ligandcoords(1)=y(3*j-2)-atomblockcentre(i,1) 
422:              ligandcoords(2)=y(3*j-1)-atomblockcentre(i,2) 
423:              ligandcoords(3)=y(3*j  )-atomblockcentre(i,3) 
424: ! rotate the block of atoms with random rotation matrix 
425:              ligandcoordsrotated=MATMUL(rotationmatrix(i,:,:),ligandcoords) 
426: ! translate back the block of atoms to the original centre of their coordinates 
427:              y(3*j-2)=ligandcoordsrotated(1)+atomblockcentre(i,1) 
428:              y(3*j-1)=ligandcoordsrotated(2)+atomblockcentre(i,2) 
429:              y(3*j  )=ligandcoordsrotated(3)+atomblockcentre(i,3) 
430:           end do             
431:        end do 
432:  
433: ! 2) RANDOM CARTESIAN PERTURBATION 
434:        if((LIGMOVET).and.(ligcartstep.gt.0)) then 
435:           do i=1,nmovableatoms 
436:             if(frozen(i)) cycle 
437: !            if(.not.enabledmove(i)) cycle 
438:             j=movableatomlist(i) 
439:             randomx=2*(DPRAND()-0.5D0) 
440:             randomy=2*(DPRAND()-0.5D0) 
441:             randomz=2*(DPRAND()-0.5D0) 
442:             y(3*j-2)=y(3*j-2)+randomx*ligcartstep*LOCALSTEP 
443:             y(3*j-1)=y(3*j-1)+randomy*ligcartstep*LOCALSTEP 
444:             y(3*j  )=y(3*j  )+randomz*ligcartstep*LOCALSTEP 
445:           enddo 
446:        endif 
447: ! 3) RIGID TRANSLATION OF THE BLOCK OF ATOMS 
448:        if((LIGMOVET).and.(ligtransstep.gt.0).and..not.BLOCKMOVET) then 
449:           randomx=2*(DPRAND()-0.5D0) 
450:           randomy=2*(DPRAND()-0.5D0) 
451:           randomz=2*(DPRAND()-0.5D0) 
452:           do i=1,nmovableatoms 
453:             if(frozen(i)) cycle 
454: !             if(.not.enabledmove(i)) cycle 
455:              j=movableatomlist(i) 
456:              y(3*j-2)=y(3*j-2)+randomx*ligtransstep*LOCALSTEP 
457:              y(3*j-1)=y(3*j-1)+randomy*ligtransstep*LOCALSTEP 
458:              y(3*j  )=y(3*j  )+randomz*ligtransstep*LOCALSTEP 
459:           enddo 
460:        endif 
461:        if((LIGMOVET).and.(ligtransstep.gt.0).and.BLOCKMOVET) then 
462:           offset1=0 
463: !          WRITE(*,*) 'nblocks ', nblocks 
464:           do i=1,nblocks 
465:             if(frozen(i)) cycle 
466: !             if(.not.enabledmove(i)) cycle 
467:              if(i>1) offset1=offset1+atomsinblock(i-1) 
468:              randomx=2*(DPRAND()-0.5D0) 
469:              randomy=2*(DPRAND()-0.5D0) 
470:              randomz=2*(DPRAND()-0.5D0) 
471:              do k=1,atomsinblock(i) 
472:                 j=movableatomlist(k+offset1) 
473:                 y(3*j-2)=y(3*j-2)+randomx*ligtransstep*LOCALSTEP 
474:                 y(3*j-1)=y(3*j-1)+randomy*ligtransstep*LOCALSTEP 
475:                 y(3*j  )=y(3*j  )+randomz*ligtransstep*LOCALSTEP 
476:              end do 
477:           end do 
478:        end if 
479:        ! compute the new distances between all atoms and 
480:        ! determine whether the overlap array has changed. 
481:        do i=1,natoms 
482: !          if(overlapt) exit 
483:           do j=i+1,natoms 
484:              distancematrix(i,j)=sqrt((y(3*i-2)-y(3*j-2))**2+(y(3*i-1)-y(3*j-1))**2+(y(3*i)-y(3*j))**2) 
485:              if(i<=nblocks.and.j<=nblocks) then ! distance between two macroions 
486:                 mindistance=macroiondist/2 
487:              else if((i<=nblocks.and.j>nblocks).or.(i>nblocks.and.j<=nblocks)) then ! distance between a macroion and a counterion 
488:                 mindistance=macroiondist/4 
489:              else if(i>nblocks.and.j>nblocks) then  ! distance between two counterions 
490:                 mindistance=1.0D0 
491:              end if 
492:              if(distancematrix(i,j)<mindistance)then 
493:                 overlaparray(i,j)=1 
494: !                overlapt=.true. 
495: !                enabledmove(j)=.true. 
496:              else 
497: !                enabledmove(j)=.false.  
498:                 overlaparray(i,j)=0  
499:              end if 
500:              if(overlaparraysave(i,j)-overlaparray(i,j)<0) then  
501:                 overlapt=.true. 
502: !                write(*,*) 'atoms close together after rotation, random perturbation and group translation: ', i,j,distancematrix(i,j) 
503: !                y(:)=ysave(:) 
504:                 nretries=nretries+1 
505:                 loopcounter=loopcounter+1 
506: !                exit  
507:              else 
508:                 overlapt=.false.  
509:              end if 
510:           end do 
511:        end do 
512:  
513:        if(overlapt) cycle 
514:  ! now we finally have no overlap, so check for percolation, and take a new step until we start off from a percolated structure 
515:        PERCT=.TRUE. 
516:        CALL PERC(y,NATOMS,PERCCUT,PERCT,.FALSE.,MYUNITNEW,.FALSE.) 
517: !      WRITE(*,*) 'perccut=', perccut 
518:        IF(.NOT.PERCT) then 
519: !          write(MYUNITNEW,*) ' macroion_moves> disconnected structure detected, undoing step' 
520:           y(:)=ysave(:) 
521:           overlapt=.true. 
522:           nretries=nretries+1  ! make sure we are not getting stuck in an infinite loop 
523:           loopcounter=loopcounter+1 
524:           cycle 
525:        ELSE 
526:           overlapt=.false. 
527:        END IF 
528:      END DO  ! while(overlapt) 
529:    END IF 
530: ! reset ligcartstep and ligtransstep 
531:    ligcartstep=cartstepsave 
532:    ligtransstep=transstepsave 
533:  
534: END SUBROUTINE MACROION_MOVES 
535:  
536: END MODULE MOVES263: END MODULE MOVES


r31713/rigidbaa.f90 2017-01-21 10:38:05.210800557 +0000 r31712/rigidbaa.f90 2017-01-21 10:38:06.790938489 +0000
 62: !     CT = cos(THETA) 62: !     CT = cos(THETA)
 63: !     ST = sin(THETA) 63: !     ST = sin(THETA)
 64: !     I3(3,3) = 3x3 identity matrix 64: !     I3(3,3) = 3x3 identity matrix
 65: !     E(3,3) = the skew-symmetric matrix obtained from a unit vector  65: !     E(3,3) = the skew-symmetric matrix obtained from a unit vector 
 66: !              parallel to P (equation (2) in the paper) 66: !              parallel to P (equation (2) in the paper)
 67: !     ESQ(3,3) = the square of E 67: !     ESQ(3,3) = the square of E
 68:       DOUBLE PRECISION :: P(3), PN(3), THETA, THETA2, THETA3, CT, ST, I3(3,3), E(3,3), ESQ(3,3) 68:       DOUBLE PRECISION :: P(3), PN(3), THETA, THETA2, THETA3, CT, ST, I3(3,3), E(3,3), ESQ(3,3)
 69:  69: 
 70: !     DEk = derivate of E with respect to the kth component of P 70: !     DEk = derivate of E with respect to the kth component of P
 71:       DOUBLE PRECISION :: DE1(3,3), DE2(3,3), DE3(3,3), RM(3,3), DRM1(3,3), DRM2(3,3), DRM3(3,3) 71:       DOUBLE PRECISION :: DE1(3,3), DE2(3,3), DE3(3,3), RM(3,3), DRM1(3,3), DRM2(3,3), DRM3(3,3)
 72:       LOGICAL, intent(in) :: GTEST 72:       LOGICAL          :: GTEST
 73:  73: 
 74: !     Set the values of the idenity matrix I3 74: !     Set the values of the idenity matrix I3
 75:       I3(:,:) = 0.D0 75:       I3(:,:) = 0.D0
 76:       I3(1,1) = 1.D0; I3(2,2) = 1.D0; I3(3,3) = 1.D0 76:       I3(1,1) = 1.D0; I3(2,2) = 1.D0; I3(3,3) = 1.D0
 77:  77: 
 78: !     Calculate the value of THETA2 as the square modulus of P 78: !     Calculate the value of THETA2 as the square modulus of P
 79:       THETA2  = DOT_PRODUCT(P,P) 79:       THETA2  = DOT_PRODUCT(P,P)
 80:  80: 
 81:       IF (THETA2 < 1.0D-12) THEN 81:       IF (THETA2 < 1.0D-12) THEN
 82: !        Execute if the angle of rotation is zero 82: !        Execute if the angle of rotation is zero


r31713/sandbox.f90 2017-01-21 10:38:05.850856532 +0000 r31712/sandbox.f90 2017-01-21 10:38:07.066962518 +0000
 43:  43: 
 44:     ! to add a new potential, add a sandbox_[potential name] subroutine, update pairwise_sandbox, and update 44:     ! to add a new potential, add a sandbox_[potential name] subroutine, update pairwise_sandbox, and update
 45:     ! sandbox_input 45:     ! sandbox_input
 46:  46: 
 47:     ! table of potentials and their integer IDs func_id: 47:     ! table of potentials and their integer IDs func_id:
 48:     ! lj 1 48:     ! lj 1
 49:     ! coulomb 2 49:     ! coulomb 2
 50:     ! dipole 3 50:     ! dipole 3
 51:     ! chiro 4 51:     ! chiro 4
 52:     ! morse 5 52:     ! morse 5
 53:     ! coulomb with distance dependent dielectric 6 
 54:  
 55:  53: 
 56: contains 54: contains
 57:     subroutine initialize_sandbox_molecule(mol, mol_index) 55:     subroutine initialize_sandbox_molecule(mol, mol_index)
 58:         implicit none 56:         implicit none
 59:         type(sandbox_molecule), target, intent(inout) :: mol 57:         type(sandbox_molecule), target, intent(inout) :: mol
 60:         integer, intent(in) :: mol_index 58:         integer, intent(in) :: mol_index
 61:  59: 
 62:         integer :: site_index 60:         integer :: site_index
 63:         type(sandbox_site), pointer :: site 61:         type(sandbox_site), pointer :: site
 64:         double precision :: dummy(3, 3) 62:         double precision :: dummy(3, 3)
178:             elseif(func_id == 4) then176:             elseif(func_id == 4) then
179:                 ! chiro177:                 ! chiro
180:                 ! sigma_0 mu^2 cos_gamma sin_gamma178:                 ! sigma_0 mu^2 cos_gamma sin_gamma
181:                 call sandbox_chiro(gtest, params(1), params(2), params(3), params(4), sitei%r - sitej%r, sitei%mu_hat, &179:                 call sandbox_chiro(gtest, params(1), params(2), params(3), params(4), sitei%r - sitej%r, sitei%mu_hat, &
182:                     & sitej%mu_hat, sitei%mu_hat_deriv, sitej%mu_hat_deriv, energy_contrib, grad_contrib)180:                     & sitej%mu_hat, sitei%mu_hat_deriv, sitej%mu_hat_deriv, energy_contrib, grad_contrib)
183:             elseif(func_id == 5) then181:             elseif(func_id == 5) then
184:                 ! morse182:                 ! morse
185:                 ! eps_0 beta r_equilibrium183:                 ! eps_0 beta r_equilibrium
186:                 call sandbox_morse(gtest, params(1), params(2), params(3), &184:                 call sandbox_morse(gtest, params(1), params(2), params(3), &
187:                     sitei%r - sitej%r, energy_contrib, grad_contrib)185:                     sitei%r - sitej%r, energy_contrib, grad_contrib)
188:             elseif(func_id == 6) then 
189:                 ! coulomb 
190:                 ! sigma_0 k_c*q*q 
191:                 call sandbox_coulomb_ddiel(gtest, params(1), params(2), sitei%r - sitej%r, energy_contrib, grad_contrib) 
192:             else186:             else
193:                 print *, "pairwise_sandbox> don't have potential with func_id='", func_id, &187:                 print *, "pairwise_sandbox> don't have potential with func_id='", func_id, &
194:                     & "' for interaction of classes", classi, classj188:                     & "' for interaction of classes", classi, classj
195:                 return189:                 return
196:             end if190:             end if
197: 191: 
198:             if(gtest) then192:             if(gtest) then
199:                 ! if site is not at molecule origin, need to compute dr^i/dp^I193:                 ! if site is not at molecule origin, need to compute dr^i/dp^I
200:                 if(any(sitei%r_molframe /= 0.D0)) then194:                 if(any(sitei%r_molframe /= 0.D0)) then
201:                     do k = 1, 3195:                     do k = 1, 3
272:         ! energy: pairwise energy266:         ! energy: pairwise energy
273:         ! gradient: (1:3) = radial for site i, (4:6) = radial j, (7:9) = angular for i, (10:12), angular j267:         ! gradient: (1:3) = radial for site i, (4:6) = radial j, (7:9) = angular for i, (10:12), angular j
274: 268: 
275:         logical, intent(in) :: gtest269:         logical, intent(in) :: gtest
276:         double precision, intent(in) :: sigma_0, kqq, rij(3)270:         double precision, intent(in) :: sigma_0, kqq, rij(3)
277:         double precision, intent(out) :: energy, grad(12)271:         double precision, intent(out) :: energy, grad(12)
278: 272: 
279:         double precision :: rijmod, rij_hat(3) ! |rij| and normal vector pointing rij273:         double precision :: rijmod, rij_hat(3) ! |rij| and normal vector pointing rij
280: 274: 
281:         ! r = (|rij| / sigma_0), r_pm1 = r to power minus 1 = r^-1, etc275:         ! r = (|rij| / sigma_0), r_pm1 = r to power minus 1 = r^-1, etc
282:         double precision :: r_pm1, r_pm2 276:         double precision :: r_pm1, r_pm2
283: 277: 
284:         energy = 0.D0278:         energy = 0.D0
285:  
286:         if(gtest) grad(:) = 0.D0279:         if(gtest) grad(:) = 0.D0
287: 280: 
288:         rijmod = sqrt(dot_product(rij, rij))281:         rijmod = sqrt(dot_product(rij, rij))
289: 282: 
290:         r_pm1 = sigma_0 / rijmod283:         r_pm1 = sigma_0 / rijmod
291: 284: 
292:         energy = kqq * r_pm1285:         energy = kqq * r_pm1
293: 286: 
294:         if(gtest) then287:         if(gtest) then
295:             rij_hat(:) = rij / rijmod288:             rij_hat(:) = rij / rijmod
296:             r_pm2 = r_pm1 * r_pm1289:             r_pm2 = r_pm1 * r_pm1
297: 290: 
298:             grad(1: 3) = -kqq / sigma_0 * r_pm2 * rij_hat(:) 291:             grad(1: 3) = -kqq / sigma_0 * r_pm2 * rij_hat(:)
299:             grad(4: 6) = -grad(1: 3)292:             grad(4: 6) = -grad(1: 3)
300:             grad(7: 9) = 0.D0293:             grad(7: 9) = 0.D0
301:             grad(10: 12) = 0.D0294:             grad(10: 12) = 0.D0
302:         end if295:         end if
303: 296: 
304:     end subroutine sandbox_coulomb297:     end subroutine sandbox_coulomb
305: 298: 
306:     subroutine sandbox_coulomb_ddiel(gtest, sigma_0, kqq, rij, energy, grad) 
307:         ! gtest: if gradients are computed 
308:         ! sigma_0: length scale for rij 
309:         ! kqq: Coulomb's constant * charge_1 * charge_2 
310:         ! rij: separation between sites 
311:  
312:         ! energy: pairwise energy 
313:         ! gradient: (1:3) = radial for site i, (4:6) = radial j, (7:9) = angular for i, (10:12), angular j 
314:  
315:         logical, intent(in) :: gtest 
316:         double precision, intent(in) :: sigma_0, kqq, rij(3) 
317:         double precision, intent(out) :: energy, grad(12) 
318:  
319:         double precision :: rijmod, rij_hat(3) ! |rij| and normal vector pointing rij 
320:         double precision :: eps, epsilon_r ! epsilon_r is a distance dependent dielectric (eps*rij) 
321:  
322:  
323:         ! r = (|rij| / sigma_0), r_pm1 = r to power minus 1 = r^-1, etc 
324:         double precision :: r_pm1, r_pm2, eps_r_pm1 
325:  
326:         energy = 0.D0 
327:  
328:         eps = 4.D0 
329:  
330:         if(gtest) grad(:) = 0.D0 
331:  
332:         rijmod = sqrt(dot_product(rij, rij)) 
333:  
334:         r_pm1 = 1.0 / rijmod 
335:  
336: ! distance dependent dielectric 
337:         epsilon_r = eps * rijmod 
338:  
339:         eps_r_pm1 = 1 / epsilon_r 
340:  
341:         energy = kqq * r_pm1 * r_pm1 * 1/eps 
342:  
343:         if(gtest) then 
344:             rij_hat(:) = rij / rijmod 
345:             r_pm2 = r_pm1 * r_pm1 
346:  
347:             grad(1: 3) = -2.0D0 * kqq * r_pm2 * rij_hat(:) * r_pm1/eps 
348:             grad(4: 6) = -grad(1: 3) 
349:             grad(7: 9) = 0.D0 
350:             grad(10: 12) = 0.D0 
351:         end if 
352:  
353:     end subroutine sandbox_coulomb_ddiel 
354:  
355:     subroutine sandbox_dipole(gtest, sigma_0, mu_p2, rij, mu_hati, mu_hatj, mu_hat_derivi, mu_hat_derivj, energy, grad)299:     subroutine sandbox_dipole(gtest, sigma_0, mu_p2, rij, mu_hati, mu_hatj, mu_hat_derivi, mu_hat_derivj, energy, grad)
356:         logical, intent(in) :: gtest300:         logical, intent(in) :: gtest
357:         double precision, intent(in) :: sigma_0, mu_p2, rij(3), mu_hati(3), mu_hatj(3), &301:         double precision, intent(in) :: sigma_0, mu_p2, rij(3), mu_hati(3), mu_hatj(3), &
358:             & mu_hat_derivi(3, 3), mu_hat_derivj(3, 3)302:             & mu_hat_derivi(3, 3), mu_hat_derivj(3, 3)
359:         double precision, intent(out) :: energy, grad(12)303:         double precision, intent(out) :: energy, grad(12)
360: 304: 
361:         double precision :: rijmod, rij_hat(3)305:         double precision :: rijmod, rij_hat(3)
362:         double precision :: r_pm1, r_pm2, r_pm3, r_pm6306:         double precision :: r_pm1, r_pm2, r_pm3, r_pm6
363:         double precision :: mu_dot_ri, mu_dot_rj, mu_dot_mu307:         double precision :: mu_dot_ri, mu_dot_rj, mu_dot_mu
364:         integer :: k308:         integer :: k
758:             interactions(classi, classj)%params(3) = cos(interactions(classi, classj)%params(3))702:             interactions(classi, classj)%params(3) = cos(interactions(classi, classj)%params(3))
759:             classes(classi)%is_aniso = .true.703:             classes(classi)%is_aniso = .true.
760:             classes(classj)%is_aniso = .true.704:             classes(classj)%is_aniso = .true.
761:             classes(classi)%has_pole = .true.705:             classes(classi)%has_pole = .true.
762:             classes(classj)%has_pole = .true.706:             classes(classj)%has_pole = .true.
763:         else if(potential_name == 'morse') then707:         else if(potential_name == 'morse') then
764:             ! morse eps_0 beta r_equilirbium -> eps beta re708:             ! morse eps_0 beta r_equilirbium -> eps beta re
765:             interactions(classi, classj)%func_id = 5709:             interactions(classi, classj)%func_id = 5
766:             allocate(interactions(classi, classj)%params(5))710:             allocate(interactions(classi, classj)%params(5))
767:             read(sand_unit, *, iostat=stat) tmp, tmp, buffer, interactions(classi, classj)%params(1: 3)711:             read(sand_unit, *, iostat=stat) tmp, tmp, buffer, interactions(classi, classj)%params(1: 3)
768:         else if(potential_name == 'coulomb_ddiel') then 
769:             ! coulomb_ddiel sigma_0 k_c q q -> sigma_0 kcqq 
770:             interactions(classi, classj)%func_id = 6 
771:             allocate(interactions(classi, classj)%params(4)) 
772:             read(sand_unit, *, iostat=stat) tmp, tmp, buffer, interactions(classi, classj)%params(1: 4) 
773:             interactions(classi, classj)%params(2) = interactions(classi, classj)%params(2) * & 
774:                 & interactions(classi, classj)%params(3) * interactions(classi, classj)%params(4) 
775:         else712:         else
776:             print *, 'sandbox_input> entered unknown interaction name!'713:             print *, 'sandbox_input> entered unknown interaction name!'
777:             exit714:             exit
778:         end if715:         end if
779:     end do716:     end do
780: 717: 
781:     ! initialize molecules718:     ! initialize molecules
782:     do mol_index = 1, size(molecules)719:     do mol_index = 1, size(molecules)
783:         call initialize_sandbox_molecule(molecules(mol_index), mol_index)720:         call initialize_sandbox_molecule(molecules(mol_index), mol_index)
784:     end do721:     end do


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0