hdiff output

r32854/gay-berne.f90 2017-06-26 22:30:10.706450026 +0100 r32853/gay-berne.f90 2017-06-26 22:30:12.070468206 +0100
5292:       MG = (1.D0 - LAMDA) * AEINV + LAMDA * BEINV5292:       MG = (1.D0 - LAMDA) * AEINV + LAMDA * BEINV
5293: 5293: 
5294:       CALL MTRXIN(MG, MGINV)5294:       CALL MTRXIN(MG, MGINV)
5295: 5295: 
5296:       MGINVR  =  MATMUL(MGINV, RIJ) 5296:       MGINVR  =  MATMUL(MGINV, RIJ) 
5297: 5297: 
5298:       SLMD =  - LAMDA * (1.D0 - LAMDA) * DOT_PRODUCT(RIJ,MGINVR)5298:       SLMD =  - LAMDA * (1.D0 - LAMDA) * DOT_PRODUCT(RIJ,MGINVR)
5299: 5299: 
5300:       RETURN5300:       RETURN
5301:       END SUBROUTINE OBJCTF5301:       END SUBROUTINE OBJCTF
5302:  
5303: SUBROUTINE GENERATECGDIMER 
5304: ! this will be useful if we have a reference dimer and would like to get the 
5305: ! coarse-grained coordinates for it (virus capsid protein dimers, for example), 
5306: ! by starting from the angle-axis information about the rigid monomer containing 
5307: ! an arbitrary number of ellipsoids. For this, we need as input the coordinate centres for both units, 
5308: ! and the rotation matrix describing the rotation of the first unit to the second one. Output  
5309: ! will be the angle-axis coordinates of the second unit. 
5310: ! Individual orientations of ellipsoids should go in the pysites.xyz file. 
5311: use rotations 
5312: implicit none 
5313: integer i,j,k, realnatoms, nlines 
5314: double precision :: y(12),atomblockcentre(2,3), rotationmatrix(3,3) 
5315:  
5316: open(255,file='coords',status='old') 
5317:  
5318: call determinelines(255,nlines) 
5319: if(nlines>4) write(*,*) 'ERROR - we should be generating a rigid body dimer, so we need 4 coordinates in the coords file' 
5320:  
5321: do i=1,nlines 
5322:    read(255,*) y(3*i-2), y(3*i-1), y(3*i) 
5323: end do 
5324:  
5325: close(255) 
5326:  
5327: open(256,file='rotationmatrix',status='old') 
5328:  do j=1,3 
5329:    read(256,*) rotationmatrix(1:3,j) 
5330:  end do 
5331:  write(*,*) 'read rotation matrix:' 
5332:  write(*,*) rotationmatrix(:,:) 
5333:  realnatoms=nlines/2 
5334:     y(4:6)=MATMUL(y(1:3),rotationmatrix)+y(4:6)  
5335:     y(3*2+3*REALNATOMS-2:3*2+3*REALNATOMS)=rot_rotate_aa(y(3*2+3*REALNATOMS-2:3*2+3*REALNATOMS),rot_mx2aa(rotationmatrix(:,:))) 
5336: write(*,*) 'finished generating dimer coordinates' 
5337: do i=1,nlines 
5338:    write(*,*) y(3*i-2), y(3*i-1), y(3*i) 
5339: end do 
5340:  
5341: contains 
5342: ! Returns the inverse of a matrix calculated by finding the LU 
5343: ! decomposition.  Depends on LAPACK. 
5344: function inv(A) result(Ainv) 
5345:   double precision, dimension(:,:), intent(in) :: A 
5346:   double precision, dimension(size(A,1),size(A,2)) :: Ainv 
5347:  
5348:   double precision, dimension(size(A,1)) :: work  ! work array for LAPACK 
5349:   integer, dimension(size(A,1)) :: ipiv   ! pivot indices 
5350:   integer :: n, info 
5351:  
5352:   ! External procedures defined in LAPACK 
5353:   external DGETRF 
5354:   external DGETRI 
5355:  
5356:   ! Store A in Ainv to prevent it from being overwritten by LAPACK 
5357:   Ainv = A 
5358:   n = size(A,1) 
5359:  
5360:   ! DGETRF computes an LU factorization of a general M-by-N matrix A 
5361:   ! using partial pivoting with row interchanges. 
5362:   call DGETRF(n, n, Ainv, n, ipiv, info) 
5363:  
5364:   if (info /= 0) then 
5365:      stop 'Matrix is numerically singular!' 
5366:   end if 
5367:  
5368:   ! DGETRI computes the inverse of a matrix using the LU factorization 
5369:   ! computed by DGETRF. 
5370:   call DGETRI(n, Ainv, n, ipiv, work, n, info) 
5371:  
5372:   if (info /= 0) then 
5373:      stop 'Matrix inversion failed!' 
5374:   end if 
5375: end function inv 
5376:  
5377:  
5378: END SUBROUTINE GENERATECGDIMER 
5379:  
5380:  
5381: SUBROUTINE DETERMINELINES(nunit,nlines) 
5382: implicit none 
5383: integer nunit, nlines, iostatus 
5384: character(len=10) check 
5385:  
5386: REWIND(nunit) 
5387:  
5388: nlines=0 
5389: do 
5390:   IF(iostatus<0) EXIT 
5391:   nlines = nlines + 1 
5392:   READ(nunit,*,iostat=iostatus) check 
5393: !  write(*,*) check,nunit 
5394: end do 
5395:   nlines = nlines - 1 
5396:   REWIND(nunit) 
5397: RETURN 
5398:  
5399:  
5400: END SUBROUTINE DETERMINELINES 
5401:  
5402:  


r32854/keywords.f 2017-06-26 22:30:10.962453438 +0100 r32853/keywords.f 2017-06-26 22:30:12.330471672 +0100
6655:             CALL READF(SIGMAP)6655:             CALL READF(SIGMAP)
6656:          ELSE6656:          ELSE
6657:             SIGMAP=0.1D06657:             SIGMAP=0.1D0
6658:          ENDIF6658:          ENDIF
6659: 6659: 
6660:          CALL ASAOOSPRINT()6660:          CALL ASAOOSPRINT()
6661: 6661: 
6662:       ELSE IF (WORD .EQ. 'SANDBOX') THEN6662:       ELSE IF (WORD .EQ. 'SANDBOX') THEN
6663: 6663: 
6664:          SANDBOXT = .TRUE.6664:          SANDBOXT = .TRUE.
6665:          RIGID = .TRUE. 
6666: 6665: 
6667:       ELSE IF (WORD .EQ. 'SILANE') THEN6666:       ELSE IF (WORD .EQ. 'SILANE') THEN
6668: 6667: 
6669:          SILANET = .TRUE.6668:          SILANET = .TRUE.
6670:          RIGID    = .TRUE.6669:          RIGID    = .TRUE.
6671:          NRBSITES = 56670:          NRBSITES = 5
6672:          ALLOCATE(SITE(NRBSITES,3))6671:          ALLOCATE(SITE(NRBSITES,3))
6673: 6672: 
6674:          CALL DEFSILANE()6673:          CALL DEFSILANE()
6675: 6674: 
7074:                 BOXLY = BOXLY*PCUTOFF7073:                 BOXLY = BOXLY*PCUTOFF
7075:                 WRITE(MYUNIT,*) "PBC Y:",BOXLY7074:                 WRITE(MYUNIT,*) "PBC Y:",BOXLY
7076:             ENDIF7075:             ENDIF
7077:             IF (SCAN(PBC,'Zz').NE.0) THEN7076:             IF (SCAN(PBC,'Zz').NE.0) THEN
7078:                 PARAMONOVPBCZ=.TRUE.7077:                 PARAMONOVPBCZ=.TRUE.
7079:                 CALL READF(BOXLZ)7078:                 CALL READF(BOXLZ)
7080:                 BOXLZ = BOXLZ*PCUTOFF7079:                 BOXLZ = BOXLZ*PCUTOFF
7081:                 WRITE(MYUNIT,*) "PBC Z:",BOXLZ7080:                 WRITE(MYUNIT,*) "PBC Z:",BOXLZ
7082:             ENDIF7081:             ENDIF
7083:          ENDIF7082:          ENDIF
7084:       ELSE IF (WORD .EQ. 'PYCGDIMER') THEN 
7085:          CALL GENERATECGDIMER 
7086:          STOP 
7087:       ELSE IF (WORD .EQ. 'PYGRAVITY') THEN7083:       ELSE IF (WORD .EQ. 'PYGRAVITY') THEN
7088:          CALL READF(PYGRAVITYC1)7084:          CALL READF(PYGRAVITYC1)
7089:          CALL READF(PYGRAVITYC2)7085:          CALL READF(PYGRAVITYC2)
7090:          EFIELDT=.TRUE.7086:          EFIELDT=.TRUE.
7091:       ELSE IF (WORD.EQ.'PYOVERLAPTHRESH') THEN7087:       ELSE IF (WORD.EQ.'PYOVERLAPTHRESH') THEN
7092:          CALL READF(PYOVERLAPTHRESH)7088:          CALL READF(PYOVERLAPTHRESH)
7093:          WRITE(MYUNIT,'(A,F8.3)') 'keyword> ellipsoids considered to overlap for an ECF value below ', PYOVERLAPTHRESH7089:          WRITE(MYUNIT,'(A,F8.3)') 'keyword> ellipsoids considered to overlap for an ECF value below ', PYOVERLAPTHRESH
7094:          IF(NITEMS.GT.2) CALL READF(PYCFTHRESH)7090:          IF(NITEMS.GT.2) CALL READF(PYCFTHRESH)
7095:          WRITE(MYUNIT,'(A,F8.3)') 'keyword> cold fusion will be diagnosed for an ECF value below ', PYCFTHRESH7091:          WRITE(MYUNIT,'(A,F8.3)') 'keyword> cold fusion will be diagnosed for an ECF value below ', PYCFTHRESH
7096: 7092: 


r32854/minpermdist.f90 2017-06-26 22:30:11.214456797 +0100 r32853/minpermdist.f90 2017-06-26 22:30:12.586475084 +0100
 76: DOUBLE PRECISION TEMPCOORDSA(DEGFREEDOMS), TEMPCOORDSB(DEGFREEDOMS) 76: DOUBLE PRECISION TEMPCOORDSA(DEGFREEDOMS), TEMPCOORDSB(DEGFREEDOMS)
 77: DOUBLE PRECISION QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES) 77: DOUBLE PRECISION QBEST(4), SITESA(3*NTSITES), SITESB(3*NTSITES)
 78: SAVE NORBIT1, NORBIT2 78: SAVE NORBIT1, NORBIT2
 79: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS) 79: INTEGER NEWPERM(NATOMS), ALLPERM(NATOMS), SAVEPERM(NATOMS)
 80: CHARACTER(LEN=5) ZSYMSAVE 80: CHARACTER(LEN=5) ZSYMSAVE
 81: COMMON /SYS/ ZSYMSAVE 81: COMMON /SYS/ ZSYMSAVE
 82:  82: 
 83: DOUBLE PRECISION :: DINV 83: DOUBLE PRECISION :: DINV
 84:  84: 
 85: DISTANCE = 0.0D0 85: DISTANCE = 0.0D0
 86: RMAT=0.0D0 
 87: WRITE(*,*) 'RMAT', RMAT 
 88: ! sf344 temporary debug 
 89: !      CALL NEWMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGID,DEBUG,RMAT) 
 90: !      RMATBEST = RMAT 
 91: !   RETURN 
 92: ! end sf344 
 93:  86: 
 94: ! jwrm2 87: ! jwrm2
 95: IF (GTHOMSONT) THEN 88: IF (GTHOMSONT) THEN
 96:    CALL GTHOMSONMINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,DISTANCE,DIST2,RMATBEST) 89:    CALL GTHOMSONMINPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,DISTANCE,DIST2,RMATBEST)
 97:    RETURN 90:    RETURN
 98: ENDIF 91: ENDIF
 99:  92: 
100: IF (RIGIDINIT) THEN 93: IF (RIGIDINIT) THEN
101:     IF(DEBUG) THEN 94:     IF(DEBUG) THEN
102:         IF(.NOT.(ANY(ABS(COORDSA(DEGFREEDOMS+1:3*NATOMS)) .GT. 1.0E-10))) THEN 95:         IF(.NOT.(ANY(ABS(COORDSA(DEGFREEDOMS+1:3*NATOMS)) .GT. 1.0E-10))) THEN


r32854/potential.f90 2017-06-26 22:30:11.814464232 +0100 r32853/potential.f90 2017-06-26 22:30:12.842478496 +0100
1028:       IF (CSMDOGUIDET) THEN ! undo guiding changes1028:       IF (CSMDOGUIDET) THEN ! undo guiding changes
1029:          CSMGP=CSMGPSAVE1029:          CSMGP=CSMGPSAVE
1030:          CSMGPINDEX=CSMGPINDEXSAVE1030:          CSMGPINDEX=CSMGPINDEXSAVE
1031:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPSAVE(1:3, 1:3, 1:2*CSMGPINDEX)1031:          PTGP(1:3, 1:3, 1:2*CSMGPINDEX)=PTGPSAVE(1:3, 1:3, 1:2*CSMGPINDEX)
1032:          CSMNORM=CSMNORMSAVE1032:          CSMNORM=CSMNORMSAVE
1033:       END IF1033:       END IF
1034: 1034: 
1035:    ELSE IF (PERMOPT .OR. PERMINVOPT .OR. DISTOPT) THEN1035:    ELSE IF (PERMOPT .OR. PERMINVOPT .OR. DISTOPT) THEN
1036: !  EREAL is the distance in this case1036: !  EREAL is the distance in this case
1037:       CALL MINPERMDIST(FIN, X, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)1037:       CALL MINPERMDIST(FIN, X, NATOMS, DEBUG, BOXLX, BOXLY, BOXLZ, PERIODIC, TWOD, EREAL, DIST2, RIGID, RMAT)
1038:       IF(DEBUG) THEN 1038: 
1039:         WRITE(MYUNIT,'(A)') 'Rotation matrix:' 
1040:         WRITE(MYUNIT,*) RMAT(:,:) 
1041:       END IF 
1042:        
1043:    ELSE IF (MKTRAPT) THEN1039:    ELSE IF (MKTRAPT) THEN
1044:       CALL MKTRAP(NATOMS, X, GRAD, EREAL)1040:       CALL MKTRAP(NATOMS, X, GRAD, EREAL)
1045:    ELSE IF (BLJCLUSTER) THEN1041:    ELSE IF (BLJCLUSTER) THEN
1046:       CALL RAD(X, GRAD, EREAL, GRADT)1042:       CALL RAD(X, GRAD, EREAL, GRADT)
1047:       CALL LJPSHIFTBINC(X, GRAD, EREAL, GRADT, SECT)1043:       CALL LJPSHIFTBINC(X, GRAD, EREAL, GRADT, SECT)
1048: 1044: 
1049:    ELSE IF (BLJCLUSTER_NOCUT) THEN1045:    ELSE IF (BLJCLUSTER_NOCUT) THEN
1050: ! ds656> Binary LJ without cutoff1046: ! ds656> Binary LJ without cutoff
1051:       CALL RAD(X, GRAD, EREAL, GRADT)1047:       CALL RAD(X, GRAD, EREAL, GRADT)
1052:       CALL BLJ_CLUST(X, GRAD, EREAL, GRADT)1048:       CALL BLJ_CLUST(X, GRAD, EREAL, GRADT)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0