hdiff output

r31547/CMakeLists.txt 2017-01-21 10:35:42.250250000 +0000 r31546/CMakeLists.txt 2017-01-21 10:35:45.534250000 +0000
  1: ########################################################### 
  2: # 
  3: # CMAKE FILE TO GENERATE libOPEP.a  
  4: #  
  5: ############################################################ 
  6:  
  7: include_directories(${CMAKE_BINARY_DIR})  1: include_directories(${CMAKE_BINARY_DIR})
  8: file(GLOB OPEP_SOURCES *.f90 *.f *.F *.F90)  2: file(GLOB OPEP_SOURCES *.f90 *.f *.F *.F90)
  9:  
 10: # Remove the interface file 
 11: file(GLOB NOT_OPEP_SOURCES opep_interface.F90) 
 12: list(REMOVE_ITEM OPEP_SOURCES ${NOT_OPEP_SOURCES} ) 
 13:  
 14: # Create the library 
 15: add_library(OPEP ${OPEP_SOURCES})  3: add_library(OPEP ${OPEP_SOURCES})
 16: add_dependencies(OPEP gminlib)  4: add_dependencies(OPEP gminlib)


r31547/commons.f90 2017-01-21 10:35:44.066250000 +0000 r31546/commons.f90 2017-01-21 10:35:47.502250000 +0000
637:       DOUBLE PRECISION, ALLOCATABLE :: REPCUT(:), NREPCUT(:), CPREPCUT(:), REPCUTFIX(:)637:       DOUBLE PRECISION, ALLOCATABLE :: REPCUT(:), NREPCUT(:), CPREPCUT(:), REPCUTFIX(:)
638:       INTEGER, ALLOCATABLE :: NREPI(:), NREPJ(:)638:       INTEGER, ALLOCATABLE :: NREPI(:), NREPJ(:)
639:       LOGICAL, ALLOCATABLE, DIMENSION(:) :: INTFROZEN  !  MXATMS639:       LOGICAL, ALLOCATABLE, DIMENSION(:) :: INTFROZEN  !  MXATMS
640:       DOUBLE PRECISION, ALLOCATABLE ::  MLPDAT(:,:)640:       DOUBLE PRECISION, ALLOCATABLE ::  MLPDAT(:,:)
641:       INTEGER, ALLOCATABLE ::  MLPOUTCOME(:)641:       INTEGER, ALLOCATABLE ::  MLPOUTCOME(:)
642:       DOUBLE PRECISION, ALLOCATABLE ::  MLQDAT(:,:)642:       DOUBLE PRECISION, ALLOCATABLE ::  MLQDAT(:,:)
643:       INTEGER, ALLOCATABLE ::  MLQOUTCOME(:)643:       INTEGER, ALLOCATABLE ::  MLQOUTCOME(:)
644: 644: 
645:       INTEGER, DIMENSION(:,:), ALLOCATABLE :: BONDS !for QCIAMBER645:       INTEGER, DIMENSION(:,:), ALLOCATABLE :: BONDS !for QCIAMBER
646: 646: 
647: !OPEP interface 
648:       LOGICAL :: OPEPT, OPEP_RNAT 
649:  
650: END MODULE COMMONS647: END MODULE COMMONS


r31547/finalio.f90 2017-01-21 10:35:44.338250000 +0000 r31546/finalio.f90 2017-01-21 10:35:47.730250000 +0000
 16: !   along with this program; if not, write to the Free Software 16: !   along with this program; if not, write to the Free Software
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19: SUBROUTINE FINALIO 19: SUBROUTINE FINALIO
 20:   USE COMMONS 20:   USE COMMONS
 21:   USE GENRIGID, ONLY : RIGIDINIT, NRIGIDBODY, NSITEPERBODY 21:   USE GENRIGID, ONLY : RIGIDINIT, NRIGIDBODY, NSITEPERBODY
 22:   USE MODAMBER 22:   USE MODAMBER
 23:   USE MODAMBER9, ONLY : COORDS1,IH,M04,AMBFINALIO_NODE 23:   USE MODAMBER9, ONLY : COORDS1,IH,M04,AMBFINALIO_NODE
 24:   USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_FINISH, AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, & 24:   USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_FINISH, AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, &
 25:        AMBER12_WRITE_XYZ 25:        AMBER12_WRITE_XYZ
 26:   USE OPEP_INTERFACE_MOD, ONLY: OPEP_FINISH, OPEP_WRITE_PDB 
 27:   USE QMODULE 26:   USE QMODULE
 28:   USE MODCHARMM 27:   USE MODCHARMM
 29:   USE AMHGLOBALS, ONLY:NMRES,IRES 28:   USE AMHGLOBALS, ONLY:NMRES,IRES
 30:   USE BGUPMOD 29:   USE BGUPMOD
 31:   USE PERMU 30:   USE PERMU
 32:   USE MODHESS, ONLY : HESS, MASSWT 31:   USE MODHESS, ONLY : HESS, MASSWT
 33:    32:   
 34:   IMPLICIT NONE 33:   IMPLICIT NONE
 35:    34:   
 36:   !   MCP 35:   !   MCP
640:            END IF639:            END IF
641:            DO J2=1,NATOMS640:            DO J2=1,NATOMS
642:               WRITE(MYUNIT3,'(3F28.20)') QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)641:               WRITE(MYUNIT3,'(3F28.20)') QMINP(J1,3*(J2-1)+1),QMINP(J1,3*(J2-1)+2),QMINP(J1,3*(J2-1)+3)
643:            ENDDO642:            ENDDO
644:            CLOSE(MYUNIT3)643:            CLOSE(MYUNIT3)
645:         ELSE644:         ELSE
646:            DO I1=1,NATOMS645:            DO I1=1,NATOMS
647:               WRITE(MYUNIT2,'(A2,3F20.10)') ih(m04+I1-1),COORDS1(3*I1-2),COORDS1(3*I1-1),COORDS1(3*I1)646:               WRITE(MYUNIT2,'(A2,3F20.10)') ih(m04+I1-1),COORDS1(3*I1-2),COORDS1(3*I1-1),COORDS1(3*I1)
648:            ENDDO647:            ENDDO
649:         ENDIF648:         ENDIF
650:       
651:      ELSE IF (OPEPT) THEN 
652:         !we only have pdb files to write 
653:         WRITE(J1_STRING,'(I6)') J1 
654:         IF (MPIT.OR.GBHT) THEN 
655:           WRITE(MYNODE_STRING,'(I6)') MYNODE + 1 
656:           CALL OPEP_WRITE_PDB(NATOMS, QMINP(J1,:),'lowest.'//TRIM(ADJUSTL(J1_STRING))//& 
657:                    &'.'//TRIM(ADJUSTL(MYNODE_STRING))//'.pdb') 
658:         ELSE 
659:           CALL OPEP_WRITE_PDB(NATOMS, QMINP(J1,:),'lowest.'//TRIM(ADJUSTL(J1_STRING))//& 
660:                    &'.pdb') 
661:         END IF 
662:         649:         
663:      ELSE IF (CHIROT) THEN650:      ELSE IF (CHIROT) THEN
664:         CALL CHIRO_OUTPUT651:         CALL CHIRO_OUTPUT
665:         652:         
666:      ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT) THEN653:      ELSE IF (ELLIPSOIDT.OR.LJCAPSIDT.OR.GBT.OR.GBDT) THEN
667:         !654:         !
668:         ! determine the centre of coordinates, then centre them655:         ! determine the centre of coordinates, then centre them
669:         !656:         !
670:         SUMX=0.0D0657:         SUMX=0.0D0
671:         SUMY=0.0d0658:         SUMY=0.0d0
1659:      CLOSE(LUNIT)1646:      CLOSE(LUNIT)
1660:      1647:      
1661:   ENDIF1648:   ENDIF
1662:   1649:   
1663:   IF(ALLOCATED(DBNAME)) DEALLOCATE(DBNAME)1650:   IF(ALLOCATED(DBNAME)) DEALLOCATE(DBNAME)
1664:   1651:   
1665:   IF (AMBER12T) THEN1652:   IF (AMBER12T) THEN
1666:      CALL AMBER12_FINISH()1653:      CALL AMBER12_FINISH()
1667:   END IF1654:   END IF
1668:   1655:   
1669:   IF (OPEPT) CALL OPEP_FINISH() 
1670:    
1671:   CALL CPU_TIME(TEND)1656:   CALL CPU_TIME(TEND)
1672:   WRITE(MYUNIT,"(A,F18.1,A)") "time elapsed ", TEND - TSTART, " seconds"1657:   WRITE(MYUNIT,"(A,F18.1,A)") "time elapsed ", TEND - TSTART, " seconds"
1673:   WRITE(MYUNIT,"(A,I18)") "Number of potential calls ", NPCALL1658:   WRITE(MYUNIT,"(A,I18)") "Number of potential calls ", NPCALL
1674:   RETURN1659:   RETURN
1675:   1660:   
1676: END SUBROUTINE FINALIO1661: END SUBROUTINE FINALIO
1677: 1662: 
1678: SUBROUTINE AMBERDUMP(J1,QMINP)1663: SUBROUTINE AMBERDUMP(J1,QMINP)
1679:     USE COMMONS1664:     USE COMMONS
1680:     USE MODAMBER1665:     USE MODAMBER


r31547/io1.f 2017-01-21 10:35:44.578250000 +0000 r31546/io1.f 2017-01-21 10:35:47.958250000 +0000
 55:          NATOMS=NATOMS+1 55:          NATOMS=NATOMS+1
 56:          READ(91,*) COORDS(3*(NATOMS-1)+1,1),COORDS(3*(NATOMS-1)+2,1),COORDS(3*(NATOMS-1)+3,1) 56:          READ(91,*) COORDS(3*(NATOMS-1)+1,1),COORDS(3*(NATOMS-1)+2,1),COORDS(3*(NATOMS-1)+3,1)
 57:          READ(91,'(A1)') DUMMY 57:          READ(91,'(A1)') DUMMY
 58:          READ(91,'(A1)') DUMMY 58:          READ(91,'(A1)') DUMMY
 59: !        WRITE(MYUNIT,'(3G20.10)') COORDS(3*(NATOMS-1)+1,1),COORDS(3*(NATOMS-1)+2,1),COORDS(3*(NATOMS-1)+3,1) 59: !        WRITE(MYUNIT,'(3G20.10)') COORDS(3*(NATOMS-1)+1,1),COORDS(3*(NATOMS-1)+2,1),COORDS(3*(NATOMS-1)+3,1)
 60:          GOTO 13 60:          GOTO 13
 61: 14       CONTINUE 61: 14       CONTINUE
 62:          CLOSE(91) 62:          CLOSE(91)
 63:       ELSEIF (AMHT) THEN 63:       ELSEIF (AMHT) THEN
 64:           write(MYUNIT,'(A)')'DUMMY    ' 64:           write(MYUNIT,'(A)')'DUMMY    '
 65:       ELSEIF (.NOT.(AMBER12T.OR.AMBERT.OR.AMBER.OR.CPMD.OR.CHRMMT.OR.DMACRYST.OR.USERPOTT.OR.OPEPT)) THEN 65:       ELSEIF (.NOT.(AMBER12T.OR.AMBERT.OR.AMBER.OR.CPMD.OR.CHRMMT.OR.DMACRYST.OR.USERPOTT)) THEN
 66: ! jwrm2> Added coordsini for GTHOMSON 66: ! jwrm2> Added coordsini for GTHOMSON
 67:          IF (GTHOMSONT) THEN 67:          IF (GTHOMSONT) THEN
 68:            CALL FILE_OPEN('coordsini', COORDS_UNIT, .FALSE.) 68:            CALL FILE_OPEN('coordsini', COORDS_UNIT, .FALSE.)
 69:          ELSE 69:          ELSE
 70:            CALL FILE_OPEN('coords', COORDS_UNIT, .FALSE.) 70:            CALL FILE_OPEN('coords', COORDS_UNIT, .FALSE.)
 71:          END IF 71:          END IF
 72: !         CLOSE(7) 72: !         CLOSE(7)
 73: !         OPEN(UNIT=7,FILE='coords',STATUS='OLD') 73: !         OPEN(UNIT=7,FILE='coords',STATUS='OLD')
 74: !js850>  this step (NATOMS=NATOMS/NPAR) is moved to earlier in the program (to 74: !js850>  this step (NATOMS=NATOMS/NPAR) is moved to earlier in the program (to
 75: !COUNTATOMS), so that the allocation is done with the correct NATOMS 75: !COUNTATOMS), so that the allocation is done with the correct NATOMS


r31547/keywords.f 2017-01-21 10:35:44.814250000 +0000 r31546/keywords.f 2017-01-21 10:35:48.194250000 +0000
 43:       USE SWMOD, ONLY: SWINIT, MWINIT 43:       USE SWMOD, ONLY: SWINIT, MWINIT
 44:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS, 44:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_SETUP, AMBER12_GET_COORDS, AMBER12_ATOMS,
 45:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA 45:      &                                  AMBER12_RESIDUES, POPULATE_ATOM_DATA
 46:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 46:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 47:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 47:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 48:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,  48:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 
 49:      &     PARSE_MLJ_PARAMS 49:      &     PARSE_MLJ_PARAMS
 50:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT 50:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT
 51:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE 51:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE
 52:       USE MOLECULAR_DYNAMICS, ONLY : MDT, MD_TSTEP, MD_GAMMA, MD_NWAIT, MD_NFREQ, MD_NSTEPS 52:       USE MOLECULAR_DYNAMICS, ONLY : MDT, MD_TSTEP, MD_GAMMA, MD_NWAIT, MD_NFREQ, MD_NSTEPS
 53:       USE OPEP_INTERFACE_MOD, ONLY : OPEP_INIT 
 54:        53:       
 55:       IMPLICIT NONE 54:       IMPLICIT NONE
 56:  55: 
 57:       DOUBLE PRECISION, ALLOCATABLE :: MLPMEAN(:), MLQMEAN(:) 56:       DOUBLE PRECISION, ALLOCATABLE :: MLPMEAN(:), MLQMEAN(:)
 58:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, LAST, IX, J1, JP, NPCOUNT, NDUMMY, INDEX, J2, J3, J4 57:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, LAST, IX, J1, JP, NPCOUNT, NDUMMY, INDEX, J2, J3, J4
 59:       INTEGER DATA_UNIT, FUNIT 58:       INTEGER DATA_UNIT, FUNIT
 60:       INTEGER MOVABLEATOMINDEX 59:       INTEGER MOVABLEATOMINDEX
 61:       LOGICAL CAT, YESNO, PERMFILE, CONFILE 60:       LOGICAL CAT, YESNO, PERMFILE, CONFILE
 62:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR, 61:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR,
 63:      &                NERROR, ECHO, LAST, CAT 62:      &                NERROR, ECHO, LAST, CAT
110:       INTEGER VECLINE, VECNUM, RESCOUNT109:       INTEGER VECLINE, VECNUM, RESCOUNT
111: 110: 
112: ! hk286 - DAMPED GROUP MOVES111: ! hk286 - DAMPED GROUP MOVES
113:       DOUBLE PRECISION GROUPATT112:       DOUBLE PRECISION GROUPATT
114: 113: 
115:       CHARACTER(LEN=120), DIMENSION(:,:), ALLOCATABLE :: KEY_WORDS114:       CHARACTER(LEN=120), DIMENSION(:,:), ALLOCATABLE :: KEY_WORDS
116:       DOUBLE PRECISION DUMMY1(NATOMS)115:       DOUBLE PRECISION DUMMY1(NATOMS)
117: 116: 
118:       INTEGER :: MAXNSETS117:       INTEGER :: MAXNSETS
119: 118: 
120:       CHARACTER(LEN=10) :: OPEP_DUMMY 
121:  
122: !      OPEN(5120, FILE = 'data')119: !      OPEN(5120, FILE = 'data')
123: !      CALL NEW_INPUT(5120, KEY_WORDS)120: !      CALL NEW_INPUT(5120, KEY_WORDS)
124: !      CLOSE(5120)121: !      CLOSE(5120)
125: !      PRINT *, KEY_WORDS122: !      PRINT *, KEY_WORDS
126: 123: 
127:       AAA=0124:       AAA=0
128:       AAB=0125:       AAB=0
129:       ABB=0126:       ABB=0
130:       PAA=0127:       PAA=0
131:       PAB=0128:       PAB=0
1164:       MLQPROB=.FALSE.1161:       MLQPROB=.FALSE.
1165:       MLQDONE=.FALSE.1162:       MLQDONE=.FALSE.
1166:       MLQNORM=.FALSE.1163:       MLQNORM=.FALSE.
1167:       MLQLAMBDA=0.0D01164:       MLQLAMBDA=0.0D0
1168:       MLQSTART=11165:       MLQSTART=1
1169: 1166: 
1170:       LJADDT=.FALSE.1167:       LJADDT=.FALSE.
1171:       LJADD2T=.FALSE.1168:       LJADD2T=.FALSE.
1172: 1169: 
1173:       DUMPMQT=.FALSE.1170:       DUMPMQT=.FALSE.
1174:  
1175: ! OPEP stuff 
1176:       OPEPT = .FALSE. 
1177:       OPEP_RNAT = .FALSE. 
1178:       1171:       
1179:       CALL FILE_OPEN('data', DATA_UNIT, .FALSE.)1172:       CALL FILE_OPEN('data', DATA_UNIT, .FALSE.)
1180:       1173:       
1181: !      OPEN (5,FILE='data',STATUS='OLD')1174: !      OPEN (5,FILE='data',STATUS='OLD')
1182: 1175: 
1183: !190   CALL INPUT(END,5)1176: !190   CALL INPUT(END,5)
1184: 190   CALL INPUT(END, DATA_UNIT)1177: 190   CALL INPUT(END, DATA_UNIT)
1185:       IF (.NOT. END) THEN1178:       IF (.NOT. END) THEN
1186:         CALL READU(WORD)1179:         CALL READU(WORD)
1187:       ENDIF1180:       ENDIF
4909: !  Specify 1D XY PBC potential4902: !  Specify 1D XY PBC potential
4910: !4903: !
4911:       ELSE IF (WORD.EQ.'ONEDPBC') THEN4904:       ELSE IF (WORD.EQ.'ONEDPBC') THEN
4912:          ONEDPBCT=.TRUE.4905:          ONEDPBCT=.TRUE.
4913:          IF (NITEMS.GT.1) CALL READI(NONEDAPBC)4906:          IF (NITEMS.GT.1) CALL READI(NONEDAPBC)
4914:          IF (MOD(NONEDAPBC,3).NE.0) THEN4907:          IF (MOD(NONEDAPBC,3).NE.0) THEN
4915:             WRITE(MYUNIT,'(A)') 'keywords> ERROR *** lattice dimension must be a multiple of three'4908:             WRITE(MYUNIT,'(A)') 'keywords> ERROR *** lattice dimension must be a multiple of three'
4916:             STOP4909:             STOP
4917:          ENDIF4910:          ENDIF
4918: 4911: 
4919:       ELSE IF (WORD.EQ.'OPEP') THEN 
4920:          OPEPT=.TRUE. 
4921:          CALL READI(NATOMS) 
4922:          WRITE(MYUNIT,'(A,I6)') 'keyword> OPEP number of atoms:',NATOMS 
4923:          CALL READA(OPEP_DUMMY) 
4924:          IF (OPEP_DUMMY.EQ.'RNA') THEN 
4925:             OPEP_RNAT = .TRUE. 
4926:             WRITE(MYUNIT,'(A)') 'keyword> RNA simulation using OPEP' 
4927:          ELSE 
4928:             WRITE(MYUNIT,'(A)') 'keyword> Protein simulation using OPEP' 
4929:          ENDIF 
4930:          IF(.NOT.ALLOCATED(COORDS1)) ALLOCATE(COORDS1(3*NATOMS)) 
4931:          IF(ALLOCATED(COORDS)) DEALLOCATE(COORDS) 
4932:          CALL OPEP_INIT(NATOMS, COORDS1(:),OPEP_RNAT) 
4933:          ALLOCATE(COORDS(3*NATOMS,NPAR)) 
4934:          DO J1=1,NPAR 
4935:            COORDS(:,J1) = COORDS1(:) 
4936:          END DO 
4937:  
4938:  
4939: ! hk2864912: ! hk286
4940:       ELSE IF (WORD.EQ.'OPTIMISEROTATION') THEN4913:       ELSE IF (WORD.EQ.'OPTIMISEROTATION') THEN
4941:          RIGIDOPTIMROTAT = .TRUE.4914:          RIGIDOPTIMROTAT = .TRUE.
4942:          CALL READF(OPTIMROTAVALUES(1))4915:          CALL READF(OPTIMROTAVALUES(1))
4943:          CALL READF(OPTIMROTAVALUES(2))4916:          CALL READF(OPTIMROTAVALUES(2))
4944:          CALL READF(OPTIMROTAVALUES(3))4917:          CALL READF(OPTIMROTAVALUES(3))
4945: 4918: 
4946: !4919: !
4947: ! NOT DOCUMENTED4920: ! NOT DOCUMENTED
4948: !4921: !


r31547/md_initialise.F90 2017-01-21 10:35:42.702250000 +0000 r31546/md_initialise.F90 2017-01-21 10:35:45.986250000 +0000
  1: module md_initialise  1: module md_initialise
  2:   use md_defs  2:   use md_defs
  3:   use md_utils  3:   use md_utils
  4:   use restart_module  4:   use restart_module
  5: contains  5: contains
  6:   ! Initialisation of the simulation  6:   ! Initialisation of the simulation
  7:   subroutine initialise(conformation,rnat)  7:   subroutine initialise(conformation)
  8:     implicit none  8:     implicit none
  9:       9:     
 10:     logical, intent(in) :: rnat 
 11:     logical :: flag 10:     logical :: flag
 12:     character(len=20) :: dummy 11:     character(len=20) :: dummy
 13:     integer :: ierror, i, i3, j3, it 12:     integer :: ierror, i, i3, j3, it
 14:     double precision :: temperature 13:     double precision :: temperature
 15:     type(t_conformations), intent(inout) :: conformation 14:     type(t_conformations), intent(inout) :: conformation
 16:     double precision :: redum 15:     double precision :: redum
 17:     double precision :: m1, m2 16:     double precision :: m1, m2
 18:  17: 
 19:     ! We first initialise the various parameters 18:     ! We first initialise the various parameters
 20:     dt=TIMESTEP*timeunit; 19:     dt=TIMESTEP*timeunit;
 21:     dt1=0.5d0*dt; 20:     dt1=0.5d0*dt;
 22:     dt2=0.5d0*dt*dt; 21:     dt2=0.5d0*dt*dt;
 23:      22:     
 24:     ! We then initialise the potential 23:     ! We then initialise the potential
 25:     RNA_simulation = rnat 
 26:     call initialise_potential() 24:     call initialise_potential()
 27:  25: 
 28:     ! We also include the inverse mass 26:     ! We also include the inverse mass
 29:     if (.not.allocated(invmass))  allocate(invmass(vecsize)) 27:     if (.not.allocated(invmass))  allocate(invmass(vecsize))
 30:     if (.not.allocated(dt1_mass)) allocate(dt1_mass(vecsize)) 28:     if (.not.allocated(dt1_mass)) allocate(dt1_mass(vecsize))
 31:     if (.not.allocated(dt2_mass)) allocate(dt2_mass(vecsize)) 29:     if (.not.allocated(dt2_mass)) allocate(dt2_mass(vecsize))
 32:  30: 
 33:  31: 
 34:     invmass = 1.0d0/mass 32:     invmass = 1.0d0/mass
 35:     dt2_mass = dt2 * invmass 33:     dt2_mass = dt2 * invmass
 82:       redu2_mass(it) = redum / beq2(it)  80:       redu2_mass(it) = redum / beq2(it) 
 83:    81:   
 84:       reduAA(it) = redu2_mass(it)/m1 82:       reduAA(it) = redu2_mass(it)/m1
 85:       reduBB(it) = redu2_mass(it)/m2 83:       reduBB(it) = redu2_mass(it)/m2
 86:          84:         
 87:       ibeq2(it) = 1.0d0 / beq2(it) 85:       ibeq2(it) = 1.0d0 / beq2(it)
 88:     end do 86:     end do
 89:        87:       
 90:     ! Generate the list of temperatures for the thermalization 88:     ! Generate the list of temperatures for the thermalization
 91:     ! Each new temperature is 1.5 as high as the previous one 89:     ! Each new temperature is 1.5 as high as the previous one
 92: !kr366> not needed with GMIN 90:     if (.not.allocated(temperatures)) allocate(temperatures(n_steps_thermalization))
 93:    ! if (.not.allocated(temperatures)) allocate(temperatures(n_steps_thermalization)) 
 94:      91:     
 95:    ! temperature = conformation%temperature 92:     temperature = conformation%temperature
 96:  93: 
 97:  94: 
 98:    ! do i=1, n_steps_thermalization 95:     do i=1, n_steps_thermalization
 99:    !    temperatures(n_steps_thermalization-i+1) = temperature 96:        temperatures(n_steps_thermalization-i+1) = temperature
100:    !    temperature = temperature / 1.5 97:        temperature = temperature / 1.5
101:    ! end do 98:     end do
102:      99:     
103:     ! we then take care of the velocities100:     ! we then take care of the velocities
104:    ! if (restart .eq. 'new') then            ! New set of velocities101:     if (restart .eq. 'new') then            ! New set of velocities
105:    !    call initialise_velocities(temperature)102:        call initialise_velocities(temperature)
106:    !    restart_simulation_time = 0103:        restart_simulation_time = 0
107:    ! else if (restart .eq. 'restart') then  ! Reuse the previous velocities and position104:     else if (restart .eq. 'restart') then  ! Reuse the previous velocities and position
108:    !    write(501,*) 'On relance'105:        write(501,*) 'On relance'
109:    !    call read_restart(restart_simulation_time,conformation)106:        call read_restart(restart_simulation_time,conformation)
110:    !    n_steps_thermalization = 0         ! Make sure that there is no thermalization107:        n_steps_thermalization = 0         ! Make sure that there is no thermalization
111:    ! else if (restart .eq. 'new_velocities') then  ! Generates a new set of velocities108:     else if (restart .eq. 'new_velocities') then  ! Generates a new set of velocities
112:    !    call read_restart(restart_simulation_time,conformation)109:        call read_restart(restart_simulation_time,conformation)
113:    !    call initialise_velocities(temperature)110:        call initialise_velocities(temperature)
114:    ! else111:     else
115:    !    stop 'Wrong choice for restart'112:        stop 'Wrong choice for restart'
116:    ! endif113:     endif
117: 114: 
118:    ! if (thermostat .eq. 'langevin') then115:     if (thermostat .eq. 'langevin') then
119:        ! Calculate the random force variance for the langevin thermostat116:        ! Calculate the random force variance for the langevin thermostat
120:    !    langevin_scalev = exp(-friction_coef*dt)117:        langevin_scalev = exp(-friction_coef*dt)
121:    ! endif118:     endif
122: 119: 
123: 120: 
124:     ! Check whether or not the filecounter file exists or not. If yes, reads it;121:     ! Check whether or not the filecounter file exists or not. If yes, reads it;
125:     ! if not, give a basic value to the counter. 122:     ! if not, give a basic value to the counter. 
126:     if (restart .ne. 'restart') then123:     if (restart .ne. 'restart') then
127:       inquire(file=COUNTER, exist=flag) 124:       inquire(file=COUNTER, exist=flag) 
128:       if (flag) then125:       if (flag) then
129:         open(unit=FCOUNTER,file=COUNTER,status='old',action='read',iostat=ierror)126:         open(unit=FCOUNTER,file=COUNTER,status='old',action='read',iostat=ierror)
130:         read(FCOUNTER,'(A13,I6)') dummy, mincounter127:         read(FCOUNTER,'(A13,I6)') dummy, mincounter
131:         close(FCOUNTER)128:         close(FCOUNTER)
139: 136: 
140:   subroutine initialise_potential()137:   subroutine initialise_potential()
141:     implicit none138:     implicit none
142:     139:     
143:     integer :: i140:     integer :: i
144:     double precision, dimension(vecsize) :: xpos    ! Positions ( (x1,y1,z1),(x2,y2,z2)...141:     double precision, dimension(vecsize) :: xpos    ! Positions ( (x1,y1,z1),(x2,y2,z2)...
145:     double precision, dimension(natoms)  :: amass   ! 1 over the masses142:     double precision, dimension(natoms)  :: amass   ! 1 over the masses
146:   143:   
147:     ! We now call the protein part to get the positions and force144:     ! We now call the protein part to get the positions and force
148:   145:   
 146:     inquire(file="conf_initiale_RNA.pdb", exist=RNA_simulation)
149:     if (RNA_simulation) then147:     if (RNA_simulation) then
150:       call initialise_RNA(PBC,BL,C_M,natoms,xpos,amass,atomic_type,force_scaling_factor, use_qbug)148:       call initialise_RNA(PBC,BL,C_M,natoms,xpos,amass,atomic_type,force_scaling_factor, use_qbug)
151:     else149:     else
152:       call initialise_protein(PBC,BL,C_M,natoms,xpos,amass,atomic_type,force_scaling_factor, use_qbug)150:       call initialise_protein(PBC,BL,C_M,natoms,xpos,amass,atomic_type,force_scaling_factor, use_qbug)
153:     endif151:     endif
154: 152: 
155:     ! We reorder the positions so that they are consistent with the rest of the program153:     ! We reorder the positions so that they are consistent with the rest of the program
156:     do i = 1, natoms154:     do i = 1, natoms
157:        x(i) = xpos(3*i-2)155:        x(i) = xpos(3*i-2)
158:        y(i) = xpos(3*i-1)156:        y(i) = xpos(3*i-1)
159:        z(i) = xpos(3*i)157:        z(i) = xpos(3*i)
160:     end do158:     end do
161:     pos = xpos159:     pos = xpos
162:    160:    
163: !    mass(1:natoms) = 1.0d0/amass161: !    mass(1:natoms) = 1.0d0/amass
164:     mass(1:3*natoms:3) = 1.0d0/amass162:     mass(1:3*natoms:3) = 1.0d0/amass
165:     mass(2:3*natoms:3) = 1.0d0/amass163:     mass(2:3*natoms:3) = 1.0d0/amass
166:     mass(3:3*natoms:3) = 1.0d0/amass164:     mass(3:3*natoms:3) = 1.0d0/amass
167:   165:   
168:     mass = mass * 2390.0 ! From amu to kcal/mol fs^2 / A^2 166:     mass = mass * 2390.0 ! From amu to kcal/mol fs^2 / A^2 
169:     return 
170:   end subroutine initialise_potential167:   end subroutine initialise_potential
171: 168: 
172: end module md_initialise169: end module md_initialise


r31547/mdmin.f90 2017-01-21 10:35:42.938250000 +0000 r31546/mdmin.f90 2017-01-21 10:35:46.214250000 +0000
 26:  contains 26:  contains
 27:  27: 
 28:     subroutine calcforce(scale,xa,fa,etot) 28:     subroutine calcforce(scale,xa,fa,etot)
 29:     use defs 29:     use defs
 30:     use calcforces 30:     use calcforces
 31:     double precision :: scale 31:     double precision :: scale
 32:     double precision, intent(out) :: etot 32:     double precision, intent(out) :: etot
 33:     double precision, dimension(vecsize) :: xpos, xfor 33:     double precision, dimension(vecsize) :: xpos, xfor
 34:     double precision, dimension(vecsize), intent(inout) :: xa 34:     double precision, dimension(vecsize), intent(inout) :: xa
 35:     double precision, dimension(vecsize), intent(out) :: fa 35:     double precision, dimension(vecsize), intent(out) :: fa
  36:     
 36:  37: 
 37:     ! Call the force field 38:     ! Call the force field
 38:     if (RNA_simulation) then 39:     if (RNA_simulation) then
 39:       call calcforce_RNA(scale,xa,fa,etot,.true., 0.0d0) 40:       call calcforce_RNA(scale,xa,fa,etot,.true., 0.0d0)
 40:     else 41:     else
 41:       call calcforce_protein(scale,xpos,xfor,etot) 42:       call calcforce_protein(scale,xpos,xfor,etot)
 42:     endif 43:     endif
 43:      44:     
 44:   end subroutine calcforce 45:   end subroutine calcforce
 45:    46:   


r31547/opep_interface.F90 2017-01-21 10:35:43.166250000 +0000 r31546/opep_interface.F90 2017-01-21 10:35:46.438250000 +0000
  1: MODULE OPEP_INTERFACE_MOD  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/OPEPinterface/opep_interface.F90' in revision 31546
  2: #ifndef DUMMY_OPEP 
  3:   USE DEFS 
  4:   USE MD_INITIALISE 
  5:   USE CALCFORCES 
  6: #endif 
  7:   IMPLICIT NONE 
  8:  
  9: !****************************************************************************** 
 10: ! Types and parameters for writing pdb ATOM coordinate records to pdb 
 11: ! files. 
 12:   type pdb_atom_data 
 13:     character (len = 6)    :: record_name 
 14:     integer                :: atom_number 
 15:     character (len = 4)    :: atom_name 
 16:     character (len = 1)    :: alt_loc_indicator 
 17:     character (len = 3)    :: residue_name 
 18:     character (len = 1)    :: chain_id 
 19:     integer                :: residue_number 
 20:     character (len = 1)    :: insertion_code 
 21:     double precision       :: x_coord 
 22:     double precision       :: y_coord 
 23:     double precision       :: z_coord 
 24:     double precision       :: occupancy 
 25:     double precision       :: b_factor 
 26:     character (len = 2)    :: element 
 27:     character (len = 2)    :: charge 
 28:   end type pdb_atom_data 
 29:  
 30: ! This defines the format to be used when writing ATOM lines for PDB 
 31: ! files.     
 32:   integer, parameter                   :: pdb_atom_data_size = 15 
 33:   type(pdb_atom_data), parameter       :: null_pdb_atom_data = & 
 34:     pdb_atom_data('ATOM  ',0,'    ',' ','   ',' ',0,'',0.d0,0.d0,0.d0,1.d0,& 
 35:                  &0.d0,'  ','  ') 
 36:   character (len=*), parameter         :: atom_string_format = & 
 37:   &'(a6,i5,1x,a4,a1,a3,1x,a1,i4,a1,3x,f8.3,f8.3,f8.3,f6.2,f6.2,10x,a2,a2)' 
 38: !****************************************************************************** 
 39:  
 40:   CHARACTER (LEN=4), ALLOCATABLE , SAVE :: RES_NAMES(:), ATOM_NAMES(:) 
 41:   INTEGER, ALLOCATABLE ,SAVE :: RES_NUMS(:) 
 42:  
 43: #ifndef DUMMY_OPEP 
 44:   TYPE (t_conformations), SAVE :: CONF 
 45: #endif 
 46:  
 47:   CONTAINS 
 48:    
 49:   SUBROUTINE OPEP_INIT(NATOM_,COORDS,RNAT) 
 50: #ifndef DUMMY_OPEP 
 51:     USE PORFUNCS 
 52:     USE COMMONS, ONLY: MYUNIT,DEBUG 
 53: #endif 
 54:     INTEGER, INTENT(IN) :: NATOM_ 
 55:     LOGICAL, INTENT(IN) :: RNAT 
 56:     DOUBLE PRECISION, INTENT(OUT) :: COORDS(3*NATOM_) 
 57: #ifndef DUMMY_OPEP 
 58:     !variables needed for creating everything for the output routine later on 
 59:     INTEGER :: PDB_FILE_UNIT, GETUNIT, L 
 60:     CHARACTER (LEN=10) :: WORD 
 61:     LOGICAL :: END_PDB 
 62:     NATOMS = NATOM_  !set the natoms variable for the OPEP interface - this is 
 63:                      !identical to GMIN, but has to be separated by name from 
 64:                      !GMIN as the OPEP internal definitions also use the name 
 65:                      !natoms 
 66:     !start initialising 
 67:     ALLOCATE(RES_NAMES(NATOMS)) 
 68:     ALLOCATE(RES_NUMS(NATOMS)) 
 69:     ALLOCATE(ATOM_NAMES(NATOMS)) 
 70:     IF (DEBUG) WRITE(MYUNIT,*) 'OPEP_init> call OPEP definitions and initialisation' 
 71:     CALL DEFINITIONS() !get all variables needed for OPEP 
 72:     CALL INITIALISE(CONF,RNAT) !initialise force field  
 73:    
 74:     COORDS=pos !pos is a variable from DEFS initiliased in INITIALISE 
 75:     !everything needed for calculations is set up now  
 76:     !before returning to the GMIN routines, initialise inofrmation we need for 
 77:     !the output at the end of the run (need residue and atom names) 
 78:     IF (DEBUG) WRITE(MYUNIT,*) 'OPEP_init> get all information needed for writing output' 
 79:     CALL FILE_OPEN("conf_initiale.pdb",PDB_FILE_UNIT,.FALSE.)  !open pdb   
 80: 30  CALL INPUT(END_PDB, PDB_FILE_UNIT)    
 81:     IF (END_PDB) THEN  !reached the end of the pdb, leave the GOTO loop 
 82:       GOTO 40 
 83:     ELSE 
 84:       CALL READU(WORD) !read the first entry, need to get all ATOM entries 
 85:     ENDIF 
 86:     IF (WORD .EQ. "ATOM") THEN 
 87:       !pdb needs to be correct format as specified below 
 88:       CALL READI(L)                  !second entry: atom number 
 89:       CALL READA(ATOM_NAMES(L))      !third entry: atom name 
 90:       CALL READA(RES_NAMES(L))       !fourth entry: residue name 
 91:       CALL READI( RES_NUMS(L))       !fifth entry: residue number 
 92:     END IF 
 93:     GOTO 30 
 94:     !close the pdb file and return to main program 
 95: 40  CLOSE(PDB_FILE_UNIT) 
 96:     IF (DEBUG) WRITE(MYUNIT,*) 'OPEP_init> finished potential initialisation' 
 97: #endif 
 98:     RETURN 
 99:   END SUBROUTINE OPEP_INIT 
100:  
101:  
102:   SUBROUTINE OPEP_ENERGY_AND_GRADIENT(NATOM_,COORD,GRAD,EREAL,GRADT) 
103: #ifndef DUMMY_OPEP 
104:     USE md_defs 
105: #endif 
106:     INTEGER, INTENT(IN) :: NATOM_ 
107:     DOUBLE PRECISION, INTENT(IN) :: COORD(3*NATOM_) 
108:     DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOM_), EREAL 
109:     LOGICAL, INTENT(IN) :: GRADT 
110: #ifndef DUMMY_OPEP 
111:     !call calcforce where we either go for RNA or protein calculation 
112:     !scale factor is set in calcforce as well, use dummy here 
113:     CALL CALCFORCE(1.0D0,COORD,GRAD,EREAL) 
114:     GRAD=-GRAD 
115: #endif 
116:     RETURN 
117:   END SUBROUTINE OPEP_ENERGY_AND_GRADIENT 
118:  
119:   SUBROUTINE OPEP_NUM_HESS(NATOM_,COORDS,DELTA,HESSIAN) 
120:     !input and output variable 
121:     INTEGER, INTENT(IN) :: NATOM_ 
122:     DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOM_), DELTA 
123:     DOUBLE PRECISION, INTENT(OUT) :: HESSIAN(3*NATOM_,3*NATOM_) 
124: #ifndef DUMMY_OPEP 
125:     !variables used internally 
126:     DOUBLE PRECISION, DIMENSION(3*NATOMS) :: COORDS_PLUS,COORDS_PLUS2,COORDS_MINUS,COORDS_MINUS2 
127:     DOUBLE PRECISION, DIMENSION(3*NATOMS) :: GRAD_PLUS,GRAD_PLUS2,GRAD_MINUS,GRAD_MINUS2 
128:     INTEGER :: I 
129:     DOUBLE PRECISION :: DUMMY_ENERGY 
130:    
131:     !use central finite difference with 4 terms (Richardson interpolation) 
132:     DO I = 1,3*NATOMS 
133:       COORDS_PLUS(:) = COORDS(:) 
134:       COORDS_PLUS(I) = COORDS(I) + DELTA 
135:       CALL OPEP_ENERGY_AND_GRADIENT(NATOMS,COORDS_PLUS,GRAD_PLUS,DUMMY_ENERGY,.TRUE.) 
136:       COORDS_PLUS2(:) = COORDS(:) 
137:       COORDS_PLUS2(I) = COORDS(I) + 2.0D0 * DELTA 
138:       CALL OPEP_ENERGY_AND_GRADIENT(NATOMS,COORDS_PLUS2,GRAD_PLUS2,DUMMY_ENERGY,.TRUE.) 
139:       COORDS_MINUS(:) = COORDS(:) 
140:       COORDS_MINUS(I) = COORDS(I) - DELTA 
141:       CALL OPEP_ENERGY_AND_GRADIENT(NATOMS,COORDS_MINUS,GRAD_MINUS,DUMMY_ENERGY,.TRUE.) 
142:       COORDS_MINUS2(:) = COORDS(:) 
143:       COORDS_MINUS2(I) = COORDS(I) - 2.0D0 * DELTA 
144:       CALL OPEP_ENERGY_AND_GRADIENT(NATOMS,COORDS_MINUS2,GRAD_MINUS2,DUMMY_ENERGY,.TRUE.) 
145:       HESSIAN(I,:) = (GRAD_MINUS2(:) - 8.0D0 * GRAD_MINUS(:) + 8.0D0 *GRAD_PLUS(:) - GRAD_PLUS2(:))/(12.0D0*DELTA) 
146:     END DO 
147: #endif 
148:     RETURN 
149:   END SUBROUTINE OPEP_NUM_HESS 
150:  
151:   SUBROUTINE OPEP_WRITE_PDB(NATOM_,COORDS,PDB_NAME) 
152: #ifndef DUMMY_OPEP 
153:     USE PORFUNCS 
154: #endif 
155:     INTEGER, INTENT(IN) :: NATOM_ 
156:     DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOM_) 
157:     CHARACTER (LEN=*), INTENT(IN) :: PDB_NAME 
158: #ifndef DUMMY_OPEP 
159:     INTEGER :: PDB_UNIT, GETUNIT 
160:     INTEGER :: CURR_ATOM 
161:     TYPE(pdb_atom_data) :: CURRENT_ATOM_DATA 
162:  
163:     !get unit for pdb file 
164:     PDB_UNIT = GETUNIT() 
165:     OPEN(UNIT=PDB_UNIT, FILE=TRIM(ADJUSTL(PDB_NAME)),ACTION="WRITE") 
166:     WRITE(PDB_UNIT,'(A)') "TITLE" 
167:     DO CURR_ATOM = 1,NATOMS 
168:       CURRENT_ATOM_DATA = null_pdb_atom_data !reset atom data for new atom 
169:       CURRENT_ATOM_DATA % atom_number = CURR_ATOM !atom number 
170:       CURRENT_ATOM_DATA % atom_name = ATOM_NAMES(CURR_ATOM) !atom name 
171:       CURRENT_ATOM_DATA % residue_number = RES_NUMS(CURR_ATOM) !res number 
172:       CURRENT_ATOM_DATA % residue_name = RES_NAMES(CURR_ATOM) !res name 
173:       CURRENT_ATOM_DATA % x_coord = COORDS(3 * CURR_ATOM -2) !x 
174:       CURRENT_ATOM_DATA % y_coord = COORDS(3 * CURR_ATOM -1) !y 
175:       CURRENT_ATOM_DATA % z_coord = COORDS(3 * CURR_ATOM) !z 
176:       WRITE(PDB_UNIT, FMT = atom_string_format) CURRENT_ATOM_DATA 
177:     ENDDO 
178:     WRITE(PDB_UNIT,'(A)') "END" 
179:     CLOSE(PDB_UNIT) 
180: #endif 
181:     RETURN 
182:   END SUBROUTINE OPEP_WRITE_PDB 
183:  
184:   SUBROUTINE OPEP_FINISH() 
185: #ifndef DUMMY_OPEP 
186:     CALL END_DEFINITIONS() 
187:     DEALLOCATE(RES_NAMES) 
188:     DEALLOCATE(RES_NUMS) 
189:     DEALLOCATE(ATOM_NAMES) 
190: #endif  
191:     RETURN 
192:   END SUBROUTINE OPEP_FINISH 
193: END MODULE OPEP_INTERFACE_MOD  


r31547/potential.f90 2017-01-21 10:35:45.042250000 +0000 r31546/potential.f90 2017-01-21 10:35:48.418250000 +0000
 28:                        TRANSFORMGRAD, TRANSFORMHESSIAN, TRANSFORMCTORIGID, & 28:                        TRANSFORMGRAD, TRANSFORMHESSIAN, TRANSFORMCTORIGID, &
 29:                        GENRIGID_DISTANCECHECK, GENRIGID_COMPRESS 29:                        GENRIGID_DISTANCECHECK, GENRIGID_COMPRESS
 30:    USE MULTIPOT, ONLY: MULTIPOT_CALL 30:    USE MULTIPOT, ONLY: MULTIPOT_CALL
 31:    USE QMODULE, ONLY: QMIN, QMINP  31:    USE QMODULE, ONLY: QMIN, QMINP 
 32:    USE PERMU, ONLY: FIN 32:    USE PERMU, ONLY: FIN
 33:    USE PORFUNCS 33:    USE PORFUNCS
 34:    USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ENERGY_AND_GRADIENT, POT_ENE_REC_C,& 34:    USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ENERGY_AND_GRADIENT, POT_ENE_REC_C,&
 35: &                                   AMBER12_NUM_HESS 35: &                                   AMBER12_NUM_HESS
 36:    USE MODHESS, ONLY: HESS, RBAANORMALMODET  36:    USE MODHESS, ONLY: HESS, RBAANORMALMODET 
 37:    USE MODAMBER9, ONLY: ATMASS1 37:    USE MODAMBER9, ONLY: ATMASS1
 38:    USE OPEP_INTERFACE_MOD, ONLY: OPEP_ENERGY_AND_GRADIENT, OPEP_NUM_HESS 
 39:    USE POLIRMOD, ONLY: POLIR 38:    USE POLIRMOD, ONLY: POLIR
 40:    USE MBPOLMOD, ONLY: MBPOL 39:    USE MBPOLMOD, ONLY: MBPOL
 41:    USE SWMOD, ONLY: SWTYPE 40:    USE SWMOD, ONLY: SWTYPE
 42:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS 41:    USE MODCUDALBFGS, ONLY: CUDA_ENEGRAD_WRAPPER, CUDA_NUMERICAL_HESS
 43:    USE GAUSS_MOD, ONLY: GFIELD 42:    USE GAUSS_MOD, ONLY: GFIELD
 44:    USE RAD_MOD, ONLY: RAD, RADR 43:    USE RAD_MOD, ONLY: RAD, RADR
 45:    USE PREC, ONLY: INT32, REAL64 44:    USE PREC, ONLY: INT32, REAL64
 46:    USE TWIST_MOD, ONLY: TWISTT, TWIST 45:    USE TWIST_MOD, ONLY: TWISTT, TWIST
 47:    IMPLICIT NONE 46:    IMPLICIT NONE
 48: ! Arguments 47: ! Arguments
560:             VNEW(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)559:             VNEW(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS)
561:             HESS(DEGFREEDOMS+1:3*NATOMS,:)=0.0D0560:             HESS(DEGFREEDOMS+1:3*NATOMS,:)=0.0D0
562:             HESS(:, DEGFREEDOMS+1:3*NATOMS)=0.0D0561:             HESS(:, DEGFREEDOMS+1:3*NATOMS)=0.0D0
563:             HESS(1:DEGFREEDOMS, 1:DEGFREEDOMS)=XRIGIDHESS(1:DEGFREEDOMS, 1:DEGFREEDOMS)562:             HESS(1:DEGFREEDOMS, 1:DEGFREEDOMS)=XRIGIDHESS(1:DEGFREEDOMS, 1:DEGFREEDOMS)
564:          END IF563:          END IF
565:       END IF564:       END IF
566: ! End of second derivative calculation 565: ! End of second derivative calculation 
567: 566: 
568: ! End of AMBERT section567: ! End of AMBERT section
569: 568: 
570: !kr366> OPEP potential 
571:    ELSE IF (OPEPT) THEN 
572:       ! If coordinates are atomistic (i.e. no rigid bodies) just call the energy 
573:       ! and gradient 
574:       IF (ATOMRIGIDCOORDT) THEN 
575:          CALL OPEP_ENERGY_AND_GRADIENT(NATOMS, X, GRADATOMS, EREAL,GRADT) 
576:          GRAD(1:3*NATOMS)=GRADATOMS(:) 
577:       ! If the coordinates include rigid bodies, transform them back to being 
578:       ! atomistic first 
579:       ELSE 
580:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS) 
581:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS) 
582:          CALL OPEP_ENERGY_AND_GRADIENT(NATOMS, X, GRADATOMS, EREAL,GRADT) 
583:          GRAD(1:3*NATOMS)=GRADATOMS(:) 
584:          !Transform the gradient and coordinates back to the rigid body representation  
585:          CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD) 
586:          X(1:DEGFREEDOMS)=XRIGIDCOORDS(1:DEGFREEDOMS) 
587:          X(DEGFREEDOMS+1:3*NATOMS)=0.0D0 
588:          GRAD(1:DEGFREEDOMS)=XRIGIDGRAD(1:DEGFREEDOMS) 
589:          GRAD(DEGFREEDOMS+1:3*NATOMS)=0.0D0 
590:       END IF 
591:       !copied from OPTIM, rbody part not tested 
592:       IF (SECT) THEN 
593:          IF (ALLOCATED(HESS)) DEALLOCATE(HESS) 
594:          ALLOCATE(HESS(3*NATOMS, 3*NATOMS)) 
595:          CALL OPEP_NUM_HESS(NATOMS,X, DELTA=1.0D-5, HESSIAN=HESS(:, :)) 
596:          IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT) ) THEN 
597:             CALL TRANSFORMHESSIAN(HESS, GRADATOMS, XRIGIDCOORDS,XRIGIDHESS,RBAANORMALMODET) 
598:             HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D0 
599:             HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D0 
600:             HESS(1:DEGFREEDOMS,1:DEGFREEDOMS) =XRIGIDHESS(1:DEGFREEDOMS,1:DEGFREEDOMS) 
601:          END IF 
602:       END IF 
603:  
604: ! hk286569: ! hk286
605:    ELSE IF (CHRMMT) THEN570:    ELSE IF (CHRMMT) THEN
606: ! hk286 > Generalised rigid body571: ! hk286 > Generalised rigid body
607: ! hk286 > If rigid coords is used, then need to convert to atom coords first, 572: ! hk286 > If rigid coords is used, then need to convert to atom coords first, 
608: ! hk286 > compute energies, then convert gradients & coords back to rigid573: ! hk286 > compute energies, then convert gradients & coords back to rigid
609:       IF (ATOMRIGIDCOORDT) THEN574:       IF (ATOMRIGIDCOORDT) THEN
610:          CALL OCHARMM(X, GRAD, EREAL, GRADT)575:          CALL OCHARMM(X, GRAD, EREAL, GRADT)
611:       ELSE576:       ELSE
612:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)577:          XRIGIDCOORDS(1:DEGFREEDOMS)=X(1:DEGFREEDOMS)
613:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)578:          CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, XCOORDS, XRIGIDCOORDS)


r31547/protein-2006.f 2017-01-21 10:35:43.394250000 +0000 r31546/protein-2006.f 2017-01-21 10:35:46.666250000 +0000
  9: C     PDB FILE: (UNIT 14).  9: C     PDB FILE: (UNIT 14).
 10: C 10: C
 11: C   2.READ THE PARAMETERS FOR THE SIDE-CHAIN SIDE-CHAIN INTERACTIONS 11: C   2.READ THE PARAMETERS FOR THE SIDE-CHAIN SIDE-CHAIN INTERACTIONS
 12: C     (UNIT 28). 12: C     (UNIT 28).
 13: C 13: C
 14: C   3.COMPUTE ENERGY and FORCES (NO CUTOFF FOR THE NONBONDED-INTERACTIONS.  14: C   3.COMPUTE ENERGY and FORCES (NO CUTOFF FOR THE NONBONDED-INTERACTIONS. 
 15: C 15: C
 16: C     THIS PROGRAM WAS WRITTEN BY P. DERREUMAUX while at IBPC, UPR 9080 16: C     THIS PROGRAM WAS WRITTEN BY P. DERREUMAUX while at IBPC, UPR 9080
 17: C     CNRS, PARIS, FRANCE.  DECEMBER 1998 17: C     CNRS, PARIS, FRANCE.  DECEMBER 1998
 18: C 18: C
 19:       USE COMMONS, ONLY: MYUNIT 19: 
 20:       use ion_pair                 !! Added by YC for OPEPv5         20:       use ion_pair                 !! Added by YC for OPEPv5        
 21:  21: 
 22:       parameter (MAXPRE = 1500)    !! maximum number of residus  22:       parameter (MAXPRE = 1500)    !! maximum number of residus 
 23:       parameter (MAXNAT = MAXPRE*6)  !! maximum number of atoms 23:       parameter (MAXNAT = MAXPRE*6)  !! maximum number of atoms
 24:       parameter (MAXXC = 3*MAXNAT)  !! maximum number of cart coord  24:       parameter (MAXXC = 3*MAXNAT)  !! maximum number of cart coord 
 25:       parameter (MAXPNB = 3*MAXPRE*MAXPRE)!! max number of SC-SC interactions  25:       parameter (MAXPNB = 3*MAXPRE*MAXPRE)!! max number of SC-SC interactions 
 26:       parameter (MAXBO  = MAXNAT)  !! maximum number of bonds 26:       parameter (MAXBO  = MAXNAT)  !! maximum number of bonds
 27:       parameter (MAXTH = MAXNAT*3)  !! maximum number of bond angles  27:       parameter (MAXTH = MAXNAT*3)  !! maximum number of bond angles 
 28:       parameter (MAXPHI = MAXNAT*4)  !! maximum number of torsional angles 28:       parameter (MAXPHI = MAXNAT*4)  !! maximum number of torsional angles
 29:       parameter (MAXTTY = 50000)       !! maximum number of residue name types  29:       parameter (MAXTTY = 50000)       !! maximum number of residue name types 
128:       real*8 box_length, inv_box_length128:       real*8 box_length, inv_box_length
129: 129: 
130:       common/pbcBL/box_length, inv_box_length130:       common/pbcBL/box_length, inv_box_length
131: 131: 
132:       character(50) dummy132:       character(50) dummy
133: C________________________________________________________________________133: C________________________________________________________________________
134:       character*30 path_for_inits134:       character*30 path_for_inits
135:       logical      single_conf_init, periodicBC,CM135:       logical      single_conf_init, periodicBC,CM
136:       periodicBC = P_B_C136:       periodicBC = P_B_C
137:       CM = CofM137:       CM = CofM
138:       write(MYUNIT,*) 'PBC and CM are : ',P_B_C, CofM138:       write(*,*) 'PBC and CM are : ',P_B_C, CofM
139: C________________________________________________________________________139: C________________________________________________________________________
140: 140: 
141:       qbug = use_qbug_141:       qbug = use_qbug_
142: 142: 
143: C --- SET THE PREFACTOR FOR SCALING THE POTENTIAL - ESSENTIAL FOR MD143: C --- SET THE PREFACTOR FOR SCALING THE POTENTIAL - ESSENTIAL FOR MD
144:       scaling_factor = a_scaling_factor !!3.2 for MD with alpha and144:       scaling_factor = a_scaling_factor !!3.2 for MD with alpha and
145: c      beta OPEP propensities 145: c      beta OPEP propensities 
146: 146: 
147: C --- SET SOME PARAMETERS147: C --- SET SOME PARAMETERS
148:       SCNB = 80.0      !! divide the 1-4 VDW interactions by 8.0148:       SCNB = 80.0      !! divide the 1-4 VDW interactions by 8.0
285:       enddo285:       enddo
286:       do ih=245,264 286:       do ih=245,264 
287:          wbeta(ih-244) = score(ih) 287:          wbeta(ih-244) = score(ih) 
288:       enddo288:       enddo
289: c --  end modif for ART and MD289: c --  end modif for ART and MD
290: 290: 
291:       291:       
292: C---------------------------------------------------------------292: C---------------------------------------------------------------
293: 293: 
294:       if (natom .ne. natom_check) then294:       if (natom .ne. natom_check) then
295:         write(MYUNIT,'(2(A,I6))') 'Natom: ',natom,' Natom_check:',natom_check295:         write(*,*) 'ERROR : Number of atom does not match input file'
296:         write(MYUNIT,*) 'ERROR : Number of atom does not match input file' 
297:         stop296:         stop
298:       endif297:       endif
299: 298: 
300:       do i=1, 3*natom_check299:       do i=1, 3*natom_check
301:         xfort(i) = x(i)300:         xfort(i) = x(i)
302:       enddo301:       enddo
303: 302: 
304:       do i=1, natom_check303:       do i=1, natom_check
305:         amasses(i) = amass(i)304:         amasses(i) = amass(i)
306:         atomic_type(i) = text3(i)305:         atomic_type(i) = text3(i)


r31547/read_parameters.f90 2017-01-21 10:35:43.842250000 +0000 r31546/read_parameters.f90 2017-01-21 10:35:47.266250000 +0000
  1: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/OPEPinterface/read_parameters.f90' in revision 31546
  2: ! 
  3: !  Reads the parameter file  
  4: ! 
  5: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  6:  
  7: SUBROUTINE definitions() 
  8:     USE COMMONS, ONLY: MYUNIT 
  9:     USE defs 
 10:     USE RANDOM 
 11:     USE md_defs 
 12:     USE ion_pair 
 13:     USE PORFUNCS 
 14:    
 15:     IMPLICIT NONE 
 16:     INTEGER :: i, PARAMS_UNIT, GETUNIT 
 17:     CHARACTER(8) :: date 
 18:     CHARACTER(10) :: time 
 19:     CHARACTER(5) :: zone 
 20:     INTEGER, DIMENSION(8) :: value 
 21:     LOGICAL :: exists_already, end 
 22:     CHARACTER WORD*16 
 23:     INTEGER :: fchain 
 24:     fchain = GETUNIT()  
 25:     use_qbug = .FALSE.             !debug force field 
 26:     force_scaling_factor = 1.0D0   !scaling of potential 
 27:     ion_pair_control = .FALSE.     !ion pair potential 
 28:     ion_pair_scaling = 1.0D0       !ion potential scaling 
 29:     PBC = .FALSE.                  !periodic boundary conditions 
 30:     BL = 1.0                       !box size 
 31:     C_M = .FALSE.                  !centre of mass for pdb output 
 32:     idum = 0                       !random number seed 
 33:     usextc = .FALSE.              !use xtc and not pdb 
 34:  
 35:     PARAMS_UNIT = GETUNIT() 
 36:     CALL FILE_OPEN('OPEP_params', PARAMS_UNIT, .FALSE.) 
 37: 10  CALL INPUT(END,PARAMS_UNIT) 
 38:     IF (.NOT. END) THEN 
 39:        CALL READU(WORD) 
 40:     ELSE 
 41:        GOTO 20 
 42:     ENDIF 
 43:     IF (WORD.EQ.'    '.OR.WORD.EQ.'NOTE'.OR.WORD.EQ.'COMMENT' & 
 44:    &                  .OR.WORD.EQ.'\\'.OR.WORD.EQ."!"         & 
 45:    &                  .OR.WORD.EQ."#") THEN  
 46:          GOTO 10 
 47:     !force field debug option 
 48:     ELSEIF (WORD .EQ. 'use_qbug') THEN 
 49:        use_qbug = .TRUE.      
 50:      
 51:     !scaling factor for potential 
 52:     ELSEIF (WORD .EQ. 'Potential_Scaling_Factor') THEN 
 53:        CALL READF(force_scaling_factor) 
 54:  
 55:     !ion pair potential 
 56:     ELSEIF (WORD .EQ. 'Ion_Pair_Potential') THEN 
 57:        ion_pair_control = .TRUE. 
 58:     ELSEIF (WORD .EQ. 'Ion_Pair_Scaling') THEN 
 59:        CALL READF(ion_pair_scaling) 
 60:  
 61:     !periodic boundary condition 
 62:     ELSEIF (WORD .EQ. 'Periodic_Boundary_Condition') THEN 
 63:        PBC = .TRUE. 
 64:        CALL READF(BL) 
 65:  
 66:     !centre of mass for pdb files for periodic boundary conditions 
 67:     ELSEIF (WORD .EQ. 'PDB_center_of_mass') THEN 
 68:        C_M = .TRUE. 
 69:  
 70:     !random number seed 
 71:     ELSEIF (WORD .EQ. 'RANDOM_SEED') THEN 
 72:        CALL READI(idum) 
 73:        !use the clock to generate a random number 
 74:        IF (idum .EQ. 0) THEN                          
 75:           CALL DATE_AND_TIME(date,time,zone,value) 
 76:           idum = -1 * MOD( (1000 * value(7) + value(8)),1024) 
 77:        ENDIF 
 78:  
 79:     ELSEIF (WORD .EQ. 'usextc') THEN 
 80:        usextc = .TRUE. 
 81:     ENDIF 
 82:     GOTO 10 
 83:  
 84:   ! We now read the number of fragments after all parameters are set 
 85: 20  INQUIRE(FILE='ichain.dat',EXIST=exists_already) 
 86:     IF (exists_already) THEN 
 87:        OPEN(UNIT=fchain,FILE='ichain.dat',STATUS='UNKNOWN',ACTION='READ') 
 88:        READ(fchain,*) nfrag 
 89:        WRITE(MYUNIT,'(A,I6)') 'read ichain for restriction',nfrag      
 90:        ALLOCATE(list_fragments(nfrag,2))    
 91:        DO i=1, nfrag 
 92:           READ(fchain,*) list_fragments(i,1),list_fragments(i,2) 
 93:        ENDDO 
 94:        CLOSE(fchain) 
 95:     ELSE 
 96:        WRITE(MYUNIT,'(A)') 'ichain.dat does not exist' 
 97:        STOP 
 98:     ENDIF 
 99:  
100:     VECSIZE = 3 * NATOMS 
101:     VECSIZE1 =  VECSIZE / NFRAG 
102:  
103:  
104: !========================================================================================================= 
105:  
106:     WRITE(MYUNIT,'(A39,I12)')   ' Number of atoms                     : ', NATOMS 
107:     WRITE(MYUNIT,'(A39,I12)')   ' Number of fragments                 : ', NFRAG 
108:     WRITE(MYUNIT,'(A39,I12)')   ' Random seed                         : ', idum 
109:     WRITE(MYUNIT,'(A39,F12.6)') ' Potential scaling factor            : ', force_scaling_factor 
110:  
111:     WRITE(MYUNIT,'(A39,L12)')   ' Periodic Boundary Condition         : ', PBC 
112:     WRITE(MYUNIT,'(A39,F12.6)') ' Box Length                          : ', BL 
113:     WRITE(MYUNIT,'(A39,L12)')   ' center of mass for pdb              : ', C_M 
114:  
115:  
116:     !CALL read_parameters_md() 
117:     ALLOCATE(pos(VECSIZE))        
118:     ALLOCATE(posref(VECSIZE))        
119:     ALLOCATE(force(VECSIZE))        
120:     ALLOCATE(atomic_type(NATOMS)) 
121:     ALLOCATE(mass(vecsize)) 
122:  
123:   ! We first set-up pointers for the x, y, z components in the position and 
124:   ! forces 
125:  
126:     x    => pos(1:3*natoms:3) 
127:     y    => pos(2:3*natoms:3) 
128:     z    => pos(3:3*natoms:3) 
129:  
130:     xref => posref(1:3*natoms:3) 
131:     yref => posref(2:3*natoms:3) 
132:     zref => posref(3:3*natoms:3) 
133:  
134:     fx   => force(1:3*natoms:3) 
135:     fy   => force(2:3*natoms:3) 
136:     fz   => force(3:3*natoms:3) 
137:  
138:    RETURN 
139: END SUBROUTINE definitions 
140:  
141: !clean-up routine for end of run 
142: SUBROUTINE end_definitions() 
143:     USE defs 
144:     USE RANDOM 
145:     USE md_defs 
146:     USE ion_pair 
147:     IF (ALLOCATED(pos)) DEALLOCATE(pos) 
148:     IF (ALLOCATED(posref)) DEALLOCATE(posref) 
149:     IF (ALLOCATED(force)) DEALLOCATE(force) 
150:     IF (ALLOCATED(atomic_type)) DEALLOCATE(atomic_type) 
151:     IF (ALLOCATED(mass)) DEALLOCATE(mass) 
152:     RETURN 
153: END SUBROUTINE end_definitions 
154:  
155:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0