hdiff output

r31568/opep_interface.F90 2017-01-21 10:35:09.862250000 +0000 r31567/opep_interface.F90 2017-01-21 10:35:10.402250000 +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/OPTIM/source/OPEPinterface/opep_interface.F90' in revision 31567
  2: #ifdef __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: #ifdef __OPEP 
 44:   TYPE (t_conformations), SAVE :: CONF 
 45: #endif 
 46:  
 47:   CONTAINS 
 48:     SUBROUTINE OPEP_GET_NATOMS(NATOM_) 
 49: #ifdef __OPEP 
 50:     USE PORFUNCS 
 51:     USE COMMONS, ONLY: DEBUG 
 52:     USE KEY, ONLY: MYUNIT 
 53: #endif 
 54:     INTEGER, INTENT(OUT) :: NATOM_ 
 55: #ifdef __OPEP 
 56:     !variables needed for creating everything for the output routine later on 
 57:     INTEGER :: PDB_FILE_UNIT, GETUNIT, L 
 58:     CHARACTER (LEN=10) :: WORD 
 59:     CHARACTER (LEN=100) :: LINE 
 60:     LOGICAL :: END_PDB 
 61:     L=0 
 62:     WRITE(*,*) ' OPEP_get_atoms> parsing conf_initiale.pdb' 
 63:     CALL FILE_OPEN("conf_initiale.pdb",PDB_FILE_UNIT,.FALSE.)  !open pdb   
 64: 30  READ(PDB_FILE_UNIT,*,ERR=40,END=40) LINE 
 65:     READ(LINE,'(A4)',ERR=30) WORD 
 66:     IF (WORD .EQ. "ATOM") THEN 
 67:        L = L + 1  
 68:        GOTO 30 
 69:     ELSE IF (WORD .EQ. "END") THEN 
 70:        GOTO 40 
 71:     ELSE 
 72:        GOTO 30 
 73:     ENDIF 
 74:     !close the pdb file and return to main program 
 75: 40  CONTINUE 
 76:     CLOSE(PDB_FILE_UNIT) 
 77:     NATOM_ = L !return number of atoms 
 78:     NATOMS = L !set internal number of atoms for OPEP potential 
 79:     WRITE(*,'(A,I6)') ' OPEP_get_atoms> Number of atoms in pdb: ',L 
 80:  
 81: #endif 
 82:     RETURN 
 83:   END SUBROUTINE OPEP_GET_NATOMS 
 84:  
 85:    
 86:   SUBROUTINE OPEP_INIT(NATOM_,COORDS,ATMASS,RNAT) 
 87: #ifdef __OPEP 
 88:     USE PORFUNCS 
 89:     USE COMMONS, ONLY: DEBUG 
 90:     USE KEY, ONLY: MYUNIT 
 91: #endif 
 92:     INTEGER, INTENT(IN) :: NATOM_ 
 93:     LOGICAL, INTENT(IN) :: RNAT 
 94:     DOUBLE PRECISION, INTENT(OUT) :: COORDS(3*NATOM_), ATMASS(NATOM_) 
 95: #ifdef __OPEP 
 96:     !variables needed for creating everything for the output routine later on 
 97:     INTEGER :: PDB_FILE_UNIT, GETUNIT, L, I 
 98:     CHARACTER (LEN=10) :: WORD 
 99:     LOGICAL :: END_PDB 
100:     NATOMS = NATOM_  !set the natoms variable for the OPEP interface - this is 
101:                      !identical to OPTIM, but has to be separated by name from 
102:                      !OPTIM as the OPEP internal definitions also use the name 
103:                      !natoms 
104:     !start initialising 
105:     ALLOCATE(RES_NAMES(NATOMS)) 
106:     ALLOCATE(RES_NUMS(NATOMS)) 
107:     ALLOCATE(ATOM_NAMES(NATOMS)) 
108:     IF (DEBUG) WRITE(*,*) ' OPEP_init> call OPEP definitions and initialisation' 
109:     CALL DEFINITIONS() !get all variables needed for OPEP 
110:     CALL INITIALISE(CONF,RNAT) !initialise force field  
111:    
112:     COORDS=pos !pos is a variable from DEFS initiliased in INITIALISE 
113:     !everything needed for calculations is set up now  
114:     !before returning to the GMIN routines, initialise inofrmation we need for 
115:     !the output at the end of the run (need residue and atom names) 
116:     DO I=1,NATOMS 
117:        ATMASS(I) = mass(3*I)/2390.0 !return the mass in amu 
118:     ENDDO 
119:     IF (DEBUG) WRITE(*,*) ' OPEP_init> get all information needed for writing output' 
120:     CALL FILE_OPEN("conf_initiale.pdb",PDB_FILE_UNIT,.FALSE.)  !open pdb   
121: 30  CALL INPUT(END_PDB, PDB_FILE_UNIT)    
122:     IF (END_PDB) THEN  !reached the end of the pdb, leave the GOTO loop 
123:       GOTO 40 
124:     ELSE 
125:       CALL READU(WORD) !read the first entry, need to get all ATOM entries 
126:     ENDIF 
127:     IF (WORD .EQ. "ATOM") THEN 
128:       !pdb needs to be correct format as specified below 
129:       CALL READI(L)                  !second entry: atom number 
130:       CALL READA(ATOM_NAMES(L))      !third entry: atom name 
131:       CALL READA(RES_NAMES(L))       !fourth entry: residue name 
132:       CALL READI( RES_NUMS(L))       !fifth entry: residue number 
133:     END IF 
134:     GOTO 30 
135:     !close the pdb file and return to main program 
136: 40  CLOSE(PDB_FILE_UNIT) 
137:     IF (DEBUG) WRITE(MYUNIT,*) 'OPEP_init> finished potential initialisation' 
138: #endif 
139:     RETURN 
140:   END SUBROUTINE OPEP_INIT 
141:  
142:  
143:   SUBROUTINE OPEP_ENERGY_AND_GRADIENT(NATOM_,COORD,GRAD,EREAL,GRADT) 
144: #ifdef __OPEP 
145:     USE md_defs 
146: #endif 
147:     INTEGER, INTENT(IN) :: NATOM_ 
148:     DOUBLE PRECISION, INTENT(IN) :: COORD(3*NATOM_) 
149:     DOUBLE PRECISION, INTENT(OUT) :: GRAD(3*NATOM_), EREAL 
150:     LOGICAL, INTENT(IN) :: GRADT 
151: #ifdef __OPEP 
152:     !call calcforce where we either go for RNA or protein calculation 
153:     !scale factor is set in calcforce as well, use dummy here 
154:     CALL CALCFORCE(1.0D0,COORD,GRAD,EREAL) 
155:     GRAD=-GRAD 
156: #endif 
157:     RETURN 
158:   END SUBROUTINE OPEP_ENERGY_AND_GRADIENT 
159:  
160:   SUBROUTINE OPEP_NUM_HESS(NATOM_,COORDS,DELTA,HESSIAN) 
161:     !input and output variable 
162:     INTEGER, INTENT(IN) :: NATOM_ 
163:     DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOM_), DELTA 
164:     DOUBLE PRECISION, INTENT(OUT) :: HESSIAN(3*NATOM_,3*NATOM_) 
165: #ifdef __OPEP 
166:     !variables used internally 
167:     DOUBLE PRECISION, DIMENSION(3*NATOMS) :: COORDS_PLUS,COORDS_PLUS2,COORDS_MINUS,COORDS_MINUS2 
168:     DOUBLE PRECISION, DIMENSION(3*NATOMS) :: GRAD_PLUS,GRAD_PLUS2,GRAD_MINUS,GRAD_MINUS2 
169:     INTEGER :: I 
170:     DOUBLE PRECISION :: DUMMY_ENERGY 
171:    
172:     !use central finite difference with 4 terms (Richardson interpolation) 
173:     DO I = 1,3*NATOMS 
174:       COORDS_PLUS(:) = COORDS(:) 
175:       COORDS_PLUS(I) = COORDS(I) + DELTA 
176:       CALL OPEP_ENERGY_AND_GRADIENT(NATOMS,COORDS_PLUS,GRAD_PLUS,DUMMY_ENERGY,.TRUE.) 
177:       COORDS_PLUS2(:) = COORDS(:) 
178:       COORDS_PLUS2(I) = COORDS(I) + 2.0D0 * DELTA 
179:       CALL OPEP_ENERGY_AND_GRADIENT(NATOMS,COORDS_PLUS2,GRAD_PLUS2,DUMMY_ENERGY,.TRUE.) 
180:       COORDS_MINUS(:) = COORDS(:) 
181:       COORDS_MINUS(I) = COORDS(I) - DELTA 
182:       CALL OPEP_ENERGY_AND_GRADIENT(NATOMS,COORDS_MINUS,GRAD_MINUS,DUMMY_ENERGY,.TRUE.) 
183:       COORDS_MINUS2(:) = COORDS(:) 
184:       COORDS_MINUS2(I) = COORDS(I) - 2.0D0 * DELTA 
185:       CALL OPEP_ENERGY_AND_GRADIENT(NATOMS,COORDS_MINUS2,GRAD_MINUS2,DUMMY_ENERGY,.TRUE.) 
186:       HESSIAN(I,:) = (GRAD_MINUS2(:) - 8.0D0 * GRAD_MINUS(:) + 8.0D0 *GRAD_PLUS(:) - GRAD_PLUS2(:))/(12.0D0*DELTA) 
187:     END DO 
188: #endif 
189:     RETURN 
190:   END SUBROUTINE OPEP_NUM_HESS 
191:  
192:   SUBROUTINE OPEP_WRITE_PDB(NATOM_,COORDS,PDB_NAME) 
193: #ifdef __OPEP 
194:     USE PORFUNCS 
195: #endif 
196:     INTEGER, INTENT(IN) :: NATOM_ 
197:     DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOM_) 
198:     CHARACTER (LEN=*), INTENT(IN) :: PDB_NAME 
199: #ifdef __OPEP 
200:     INTEGER :: PDB_UNIT, GETUNIT 
201:     INTEGER :: CURR_ATOM 
202:     TYPE(pdb_atom_data) :: CURRENT_ATOM_DATA 
203:  
204:     !get unit for pdb file 
205:     PDB_UNIT = GETUNIT() 
206:     OPEN(UNIT=PDB_UNIT, FILE=TRIM(ADJUSTL(PDB_NAME)),ACTION="WRITE") 
207:     WRITE(PDB_UNIT,'(A)') "TITLE" 
208:     DO CURR_ATOM = 1,NATOMS 
209:       CURRENT_ATOM_DATA = null_pdb_atom_data !reset atom data for new atom 
210:       CURRENT_ATOM_DATA % atom_number = CURR_ATOM !atom number 
211:       CURRENT_ATOM_DATA % atom_name = ATOM_NAMES(CURR_ATOM) !atom name 
212:       CURRENT_ATOM_DATA % residue_number = RES_NUMS(CURR_ATOM) !res number 
213:       CURRENT_ATOM_DATA % residue_name = RES_NAMES(CURR_ATOM) !res name 
214:       CURRENT_ATOM_DATA % x_coord = COORDS(3 * CURR_ATOM -2) !x 
215:       CURRENT_ATOM_DATA % y_coord = COORDS(3 * CURR_ATOM -1) !y 
216:       CURRENT_ATOM_DATA % z_coord = COORDS(3 * CURR_ATOM) !z 
217:       WRITE(PDB_UNIT, FMT = atom_string_format) CURRENT_ATOM_DATA 
218:     ENDDO 
219:     WRITE(PDB_UNIT,'(A)') "END" 
220:     CLOSE(PDB_UNIT) 
221: #endif 
222:     RETURN 
223:   END SUBROUTINE OPEP_WRITE_PDB 
224:  
225:   SUBROUTINE OPEP_FINISH() 
226: #ifdef __OPEP 
227:     CALL END_DEFINITIONS() 
228:     DEALLOCATE(RES_NAMES) 
229:     DEALLOCATE(RES_NUMS) 
230:     DEALLOCATE(ATOM_NAMES) 
231: #endif  
232:     RETURN 
233:   END SUBROUTINE OPEP_FINISH 
234:  
235: END MODULE OPEP_INTERFACE_MOD  


r31568/read_parameters.f90 2017-01-21 10:35:10.150250000 +0000 r31567/read_parameters.f90 2017-01-21 10:35:10.650250000 +0000
  1: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/OPTIM/source/OPEPinterface/read_parameters.f90' in revision 31567
  2: ! 
  3: !  Reads the parameter file  
  4: ! 
  5: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
  6:  
  7: SUBROUTINE definitions() 
  8:     USE KEY, 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