hdiff output

r33360/amber_mutations.F90 2017-09-29 15:30:09.742178861 +0100 r33359/amber_mutations.F90 2017-09-29 15:30:10.414187632 +0100
116:         ENTRIES(:)=''116:         ENTRIES(:)=''
117:         CALL READ_LINE(ENTRY_,MUTATION_INFO(J1)%NENTRIES,ENTRIES)117:         CALL READ_LINE(ENTRY_,MUTATION_INFO(J1)%NENTRIES,ENTRIES)
118:         DO J2=1,MUTATION_INFO(J1)%NENTRIES118:         DO J2=1,MUTATION_INFO(J1)%NENTRIES
119:            READ(ENTRIES(J2),*) MUTATION_INFO(J1)%PROBABILITIES(J2)119:            READ(ENTRIES(J2),*) MUTATION_INFO(J1)%PROBABILITIES(J2)
120:         ENDDO120:         ENDDO
121:      ENDDO121:      ENDDO
122: 101  CONTINUE122: 101  CONTINUE
123:      CLOSE(MIUNIT)123:      CLOSE(MIUNIT)
124:      !call the grouprotation set up here (not in keywords)124:      !call the grouprotation set up here (not in keywords)
125:      CALL MUT_SETUP_GROUPROTATION(1,.FALSE.,.FALSE.,0)125:      CALL MUT_SETUP_GROUPROTATION(1,.FALSE.,.FALSE.,0)
126: !old implementation --> use file to specify residues to be rigidified 126:      IF (AMBERMUTRIGIDT) THEN
127: !     IF (AMBERMUTRIGIDT) THEN127:         CALL CREATE_RIGID_FILES()
128: !        CALL CREATE_RIGID_FILES()128:      END IF    
129: !     END IF     
130:      RETURN129:      RETURN
131:   END SUBROUTINE AMBERMUTATION_SETUP130:   END SUBROUTINE AMBERMUTATION_SETUP
132: 131: 
133:   !setup rigid bodies for AMBERMUTATIONS132:   !setup rigid bodies for AMBERMUTATIONS
134:   SUBROUTINE AMBERMUTRB_SETUP()133:   SUBROUTINE AMBERMUTRB_SETUP()
135:      IMPLICIT NONE134:      IMPLICIT NONE
136:      INTEGER :: GETUNIT, RBUNIT, J1135:      INTEGER :: GETUNIT, RBUNIT, J1
137:      LOGICAL :: YESNO136:      LOGICAL :: YESNO
138:      CHARACTER(200) :: ENTRY_137:      CHARACTER(200) :: ENTRY_
139:      CHARACTER(25) :: ENTRIES(2)138:      CHARACTER(25) :: ENTRIES(2)
159:      ENDDO158:      ENDDO
160: 159  CONTINUE159: 159  CONTINUE
161:      CLOSE(RBUNIT)160:      CLOSE(RBUNIT)
162:      RETURN161:      RETURN
163:   END SUBROUTINE AMBERMUTRB_SETUP162:   END SUBROUTINE AMBERMUTRB_SETUP
164:   163:   
165:   !mutate protein164:   !mutate protein
166:   SUBROUTINE AMBERMUT_STEP(COORDINATES , RESNUMBER)165:   SUBROUTINE AMBERMUT_STEP(COORDINATES , RESNUMBER)
167:      INTEGER , INTENT(OUT) :: RESNUMBER166:      INTEGER , INTENT(OUT) :: RESNUMBER
168:      CHARACTER(LEN=4) :: OLDRES , OLDRES1 , NEWRES , NEWRES1167:      CHARACTER(LEN=4) :: OLDRES , OLDRES1 , NEWRES , NEWRES1
169:      CHARACTER(LEN=6) :: NMUT_STRING , STARTINDEX_STRING, START_STR, SHIFT_STR168:      CHARACTER(LEN=6) :: NMUT_STRING , STARTINDEX_STRING
170:      CHARACTER(LEN=50) :: OPTION_STRING, PATH_STRING169:      CHARACTER(LEN=25) :: OPTION_STRING
171:      DOUBLE PRECISION :: COORDINATES(3*NATOMS)170:      DOUBLE PRECISION :: COORDINATES(3*NATOMS)
172:      INTEGER ::  SHIFT, START 
173:  171:  
174:      !let's store all information first in case we have to go back!172:      !let's store all information first in case we have to go back!
175:      PREVIOUS_MUTATION = MUTATION_INFO173:      PREVIOUS_MUTATION = MUTATION_INFO
176:      !we have a new mutation174:      !we have a new mutation
177:      NMUTATION = NMUTATION + 1175:      NMUTATION = NMUTATION + 1
178:      WRITE(NMUT_STRING,'(I6)') NMUTATION - 1 176:      WRITE(NMUT_STRING,'(I6)') NMUTATION - 1 
179:      !before we do anything, we save the old lowest minima177:      !before we do anything, we save the old lowest minima
180:      CALL AMBERMUT_CURR_LOWEST()178:      CALL AMBERMUT_CURR_LOWEST()
181:      !select a residue to mutate179:      !select a residue to mutate
182:      CALL SELECT_MUTATION(RESNUMBER , OLDRES1 , NEWRES1)180:      CALL SELECT_MUTATION(RESNUMBER , OLDRES1 , NEWRES1)
205: #ifdef _SVN_ROOT_203: #ifdef _SVN_ROOT_
206:      CALL SYSTEM('python '//_SVN_ROOT_//'/SCRIPTS/AMBER/BHmutation_steps/mutate_aa.py '//OLDRES//' '//NEWRES)204:      CALL SYSTEM('python '//_SVN_ROOT_//'/SCRIPTS/AMBER/BHmutation_steps/mutate_aa.py '//OLDRES//' '//NEWRES)
207:      CALL SYSTEM('python '//_SVN_ROOT_//'/SCRIPTS/AMBER/BHmutation_steps/perm_allow.py '//OPTION_STRING)205:      CALL SYSTEM('python '//_SVN_ROOT_//'/SCRIPTS/AMBER/BHmutation_steps/perm_allow.py '//OPTION_STRING)
208: #else206: #else
209:      CALL SYSTEM('python ' // mutation_script // OLDRES // ' ' // NEWRES )207:      CALL SYSTEM('python ' // mutation_script // OLDRES // ' ' // NEWRES )
210:      CALL SYSTEM('python ' // perm_allow_script //OPTION_STRING)208:      CALL SYSTEM('python ' // perm_allow_script //OPTION_STRING)
211: #endif209: #endif
212:      CALL SYSTEM('mv perm.allow perm.allow.'//TRIM(ADJUSTL(NMUT_STRING)))210:      CALL SYSTEM('mv perm.allow perm.allow.'//TRIM(ADJUSTL(NMUT_STRING)))
213:      CALL SYSTEM('mv perm.allow.new perm.allow')211:      CALL SYSTEM('mv perm.allow.new perm.allow')
214:      !create a new topology, update the residue information and adjust coordinates for unchanged residues212:      !create a new topology, update the residue information and adjust coordinates for unchanged residues
215:      CALL CREATE_NEW_TOPOLOGY(RESNUMBER ,  NEWRES , COORDS, START, SHIFT)213:      CALL CREATE_NEW_TOPOLOGY(RESNUMBER ,  NEWRES , COORDS)
216:      WRITE(START_STR,'(I6)') START 
217:      WRITE(SHIFT_STR,'(I6)') SHIFT 
218:      !create new atom groups214:      !create new atom groups
219:      IF (.NOT.AMBERMUTRIGIDT) THEN 
220: #ifdef _SVN_ROOT_215: #ifdef _SVN_ROOT_
221:         CALL SYSTEM('python ' // _SVN_ROOT_ // '/SCRIPTS/AMBER/BHmutation_steps/grouprotations.py tmp.pdb')216:      CALL SYSTEM('python ' // _SVN_ROOT_ // '/SCRIPTS/AMBER/BHmutation_steps/grouprotations.py tmp.pdb')
222: #else217: #else
223:         CALL SYSTEM('python ' // grouprotation_script // ' tmp.pdb')218:      CALL SYSTEM('python ' // grouprotation_script // ' tmp.pdb')
224: #endif219: #endif
225:      ELSE220:      IF (AMBERMUTRIGIDT) THEN
226:         OPTION_STRING='rbodyconfig '//SHIFT_STR//' '//START_STR//' tmp.pdb' 
227:         PATH_STRING='/SCRIPTS/AMBER/BHmutation_steps/rbodyconfig_grouprotations.py' 
228: #ifdef _SVN_ROOT_ 
229:          
230:         CALL SYSTEM('python '// _SVN_ROOT_ //PATH_STRING//OPTION_STRING) 
231:         CALL SYSTEM('mv rbodyconfig rbodyconfig.'//TRIM(ADJUSTL(NMUT_STRING)))221:         CALL SYSTEM('mv rbodyconfig rbodyconfig.'//TRIM(ADJUSTL(NMUT_STRING)))
232: #else 
233:         CALL SYSTEM('python ' // grouprotation_script //OPTION_STRING)       
234: #endif 
235:         CALL SYSTEM('mv rbodyconfig.new rbodyconfig.') 
236:         CALL SYSTEM('mv coordsinirigid coordsinirigid.'//TRIM(ADJUSTL(NMUT_STRING)))222:         CALL SYSTEM('mv coordsinirigid coordsinirigid.'//TRIM(ADJUSTL(NMUT_STRING)))
237:         CALL SYSTEM('cp start coordsinirigid') 
238:      ENDIF223:      ENDIF
239:      CALL SYSTEM('rm tmp.pdb')224:      CALL SYSTEM('rm tmp.pdb')
240:      !finally reinitialise AMBER with new groups, coordinates and topology225:      !finally reinitialise AMBER with new groups, coordinates and topology
241:      CALL REINITIALISE_AMBER()226:      CALL REINITIALISE_AMBER()
242:      !now remove old chiral states used for checking (the rest is done when we initialise the chirality in mc.F)227:      !now remove old chiral states used for checking (the rest is done when we initialise the chirality in mc.F)
243:      CALL DEALLOC_STATES_MUTATION()228:      CALL DEALLOC_STATES_MUTATION()
244:      RETURN229:      RETURN
245:   END SUBROUTINE AMBERMUT_STEP230:   END SUBROUTINE AMBERMUT_STEP
246: 231: 
247:   SUBROUTINE SELECT_MUTATION(RESNUMBER , OLDRES , NEWRES)232:   SUBROUTINE SELECT_MUTATION(RESNUMBER , OLDRES , NEWRES)
356:      STARTATOM = AMBER12_RESSTART(RESNUMBER)341:      STARTATOM = AMBER12_RESSTART(RESNUMBER)
357:      NRESATOM = AMBER12_RESNATOM(RESNUMBER)342:      NRESATOM = AMBER12_RESNATOM(RESNUMBER)
358:      CUNIT = GETUNIT()343:      CUNIT = GETUNIT()
359:      OPEN(UNIT=CUNIT , FILE='coords.oldres' , STATUS='NEW')344:      OPEN(UNIT=CUNIT , FILE='coords.oldres' , STATUS='NEW')
360:      DO J1 = 1,NRESATOM345:      DO J1 = 1,NRESATOM
361:         WRITE(CUNIT,'(3F20.10)') COORD(3*(STARTATOM+J1-1)-2),COORD(3*(STARTATOM+J1-1)-1),COORD(3*(STARTATOM+J1-1))346:         WRITE(CUNIT,'(3F20.10)') COORD(3*(STARTATOM+J1-1)-2),COORD(3*(STARTATOM+J1-1)-1),COORD(3*(STARTATOM+J1-1))
362:      ENDDO347:      ENDDO
363:      CLOSE(CUNIT)348:      CLOSE(CUNIT)
364:   END SUBROUTINE DUMP_RESIDUE_COORDS349:   END SUBROUTINE DUMP_RESIDUE_COORDS
365: 350: 
366:   SUBROUTINE CREATE_NEW_TOPOLOGY(RESNUMBER , NEWRES , COORDS_OLD, STARTATOM, SHIFT)351:   SUBROUTINE CREATE_NEW_TOPOLOGY(RESNUMBER , NEWRES , COORDS_OLD)
367:      INTEGER , INTENT(IN) :: RESNUMBER352:      INTEGER , INTENT(IN) :: RESNUMBER
368:      CHARACTER(LEN=4) , INTENT(IN) :: NEWRES353:      CHARACTER(LEN=4) , INTENT(IN) :: NEWRES
369:      DOUBLE PRECISION , INTENT(IN) :: COORDS_OLD(3*NATOMS)354:      DOUBLE PRECISION , INTENT(IN) :: COORDS_OLD(3*NATOMS)
370:      INTEGER, INTENT(OUT) :: STARTATOM , SHIFT 
371:      DOUBLE PRECISION , ALLOCATABLE :: COORDS_NEW(:) , COORDS_NEWRES(:,:)355:      DOUBLE PRECISION , ALLOCATABLE :: COORDS_NEW(:) , COORDS_NEWRES(:,:)
372:      DOUBLE PRECISION , ALLOCATABLE :: COORDS_RES(:)356:      DOUBLE PRECISION , ALLOCATABLE :: COORDS_RES(:)
373:      INTEGER :: J1 , TUNIT , CUNIT , CUNIT2 , GETUNIT ,  FINALATOM_OLD , FINALATOM_NEW 357:      INTEGER :: J1 , TUNIT , CUNIT , CUNIT2 , GETUNIT , STARTATOM , FINALATOM_OLD , FINALATOM_NEW , SHIFT 
374:      CHARACTER(LEN=4) :: RESNAMES(NRESIDUES)358:      CHARACTER(LEN=4) :: RESNAMES(NRESIDUES)
375: 359: 
376:      TUNIT = GETUNIT()360:      TUNIT = GETUNIT()
377:      DO J1=1,NRESIDUES361:      DO J1=1,NRESIDUES
378:         RESNAMES(J1) = AMBER12_RESNAME(J1)362:         RESNAMES(J1) = AMBER12_RESNAME(J1)
379:      ENDDO363:      ENDDO
380:      !create a leap.in file364:      !create a leap.in file
381:      OPEN(TUNIT , FILE='leap.in' , STATUS='NEW')365:      OPEN(TUNIT , FILE='leap.in' , STATUS='NEW')
382:      !currently we either go for ff14SB or ff99SB366:      !currently we either go for ff14SB or ff99SB
383:      IF (AMBERMUTFF.EQ.14) THEN367:      IF (AMBERMUTFF.EQ.14) THEN
504:         ALLOCATE(COORDS1(3*NATOMS))488:         ALLOCATE(COORDS1(3*NATOMS))
505:         IF(ALLOCATED(COORDS)) DEALLOCATE(COORDS)489:         IF(ALLOCATED(COORDS)) DEALLOCATE(COORDS)
506:         ! Read the coords from AMBER12 into COORDS1(:)490:         ! Read the coords from AMBER12 into COORDS1(:)
507:         CALL AMBER12_GET_COORDS(NATOMS, COORDS1(:))491:         CALL AMBER12_GET_COORDS(NATOMS, COORDS1(:))
508:         ALLOCATE(COORDS(3*NATOMS,NPAR))492:         ALLOCATE(COORDS(3*NATOMS,NPAR))
509:         DO J1=1,NPAR493:         DO J1=1,NPAR
510:            COORDS(:,J1) = COORDS1(:)494:            COORDS(:,J1) = COORDS1(:)
511:         END DO495:         END DO
512:         !setup the new group rotation information496:         !setup the new group rotation information
513:         CALL MUT_SETUP_GROUPROTATION(1,.FALSE.,.FALSE.,0)497:         CALL MUT_SETUP_GROUPROTATION(1,.FALSE.,.FALSE.,0)
514:         !call reinitialisation for rigid bodies498:         !setup new rigid body files and call reinitialisation
515:         IF (AMBERMUTRIGIDT) THEN499:         IF (AMBERMUTRIGIDT) THEN
516:            CALL DEALLOCATE_GENRIGID()500:            CALL DEALLOCATE_GENRIGID()
 501:            CALL CREATE_RIGID_FILES()
517:            CALL GENRIGID_READ_FROM_FILE()502:            CALL GENRIGID_READ_FROM_FILE()
518:         ENDIF503:         ENDIF
519:         !deallocate, reallocate and initialise a bunch of globals that we need to reset504:         !deallocate, reallocate and initialise a bunch of globals that we need to reset
520:         DEALLOCATE(QMINP)505:         DEALLOCATE(QMINP)
521:         ALLOCATE(QMINP(NSAVE,3*NATOMS))506:         ALLOCATE(QMINP(NSAVE,3*NATOMS))
522:         DEALLOCATE(QMINT)507:         DEALLOCATE(QMINT)
523:         ALLOCATE(QMINT(NSAVE,NATOMS))508:         ALLOCATE(QMINT(NSAVE,NATOMS))
524:         DEALLOCATE(COORDSO)509:         DEALLOCATE(COORDSO)
525:         ALLOCATE(COORDSO(3*NATOMS,NPAR))510:         ALLOCATE(COORDSO(3*NATOMS,NPAR))
526:         DEALLOCATE(VT)511:         DEALLOCATE(VT)
564:      FINALATOM_NEW = AMBER12_RESEND(RESNUMBER)549:      FINALATOM_NEW = AMBER12_RESEND(RESNUMBER)
565:      SHIFT = FINALATOM_NEW - FINALATOM_OLD550:      SHIFT = FINALATOM_NEW - FINALATOM_OLD
566:      !correct wrong information (we havent reinitialised yet, so NATOMS is still wrong)551:      !correct wrong information (we havent reinitialised yet, so NATOMS is still wrong)
567:      AMBER12_RESEND(NRESIDUES) = AMBER12_RESEND(NRESIDUES) + SHIFT552:      AMBER12_RESEND(NRESIDUES) = AMBER12_RESEND(NRESIDUES) + SHIFT
568:      AMBER12_RESNATOM(NRESIDUES) = AMBER12_RESNATOM(NRESIDUES) + SHIFT553:      AMBER12_RESNATOM(NRESIDUES) = AMBER12_RESNATOM(NRESIDUES) + SHIFT
569:      !create final input files needed554:      !create final input files needed
570:      CALL SYSTEM('cp coords.'//TRIM(ADJUSTL(NMUT_STRING))//'.1.rst coords.inpcrd')555:      CALL SYSTEM('cp coords.'//TRIM(ADJUSTL(NMUT_STRING))//'.1.rst coords.inpcrd')
571:      CALL SYSTEM('cp start.'//TRIM(ADJUSTL(NMUT_STRING))//' start')  !this one is the wrong file?556:      CALL SYSTEM('cp start.'//TRIM(ADJUSTL(NMUT_STRING))//' start')  !this one is the wrong file?
572:      CALL SYSTEM('cp atomgroups.'//TRIM(ADJUSTL(NMUT_STRING))//' atomgroups')557:      CALL SYSTEM('cp atomgroups.'//TRIM(ADJUSTL(NMUT_STRING))//' atomgroups')
573:      CALL SYSTEM('cp perm.allow.'//TRIM(ADJUSTL(NMUT_STRING))//' perm.allow')558:      CALL SYSTEM('cp perm.allow.'//TRIM(ADJUSTL(NMUT_STRING))//' perm.allow')
574:      IF (AMBERMUTRIGIDT) THEN 
575:         CALL SYSTEM('cp rbodyconfig.'//TRIM(ADJUSTL(NMUT_STRING))//' rbodyconfig')  
576:         CALL SYSTEM('cp start coordsinirigid')      
577:      ENDIF 
578:      !now reinitialise once more559:      !now reinitialise once more
579:      CALL REINITIALISE_AMBER()560:      CALL REINITIALISE_AMBER()
580:      !reset chirality561:      !reset chirality
581:      CALL DEALLOC_STATES_MUTATION()562:      CALL DEALLOC_STATES_MUTATION()
582:   END SUBROUTINE REVERSE_MUTATION563:   END SUBROUTINE REVERSE_MUTATION
583: 564: 
584:   !scoring function to judge how good mutation is565:   !scoring function to judge how good mutation is
585:   SUBROUTINE MUTATION_E(SCORE,COORDS,MODE,TERMID)566:   SUBROUTINE MUTATION_E(SCORE,COORDS,MODE,TERMID)
586:      DOUBLE PRECISION, INTENT(OUT) :: SCORE567:      DOUBLE PRECISION, INTENT(OUT) :: SCORE
587:      DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)568:      DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)


r33360/keywords.f 2017-09-29 15:30:09.970181843 +0100 r33359/keywords.f 2017-09-29 15:30:10.654190761 +0100
 27:       USE TWIST_MOD 27:       USE TWIST_MOD
 28: !       sf344> AMBER additions 28: !       sf344> AMBER additions
 29:       USE modamber9, only : coords1,amberstr,amberstr1,mdstept,inpcrd,amberenergiest, nocistransdna, nocistransrna, 29:       USE modamber9, only : coords1,amberstr,amberstr1,mdstept,inpcrd,amberenergiest, nocistransdna, nocistransrna,
 30:      &                      uachiral, ligrotscale, setchiral, STEEREDMINT, SMINATOMA, SMINATOMB, SMINK, SMINKINC, 30:      &                      uachiral, ligrotscale, setchiral, STEEREDMINT, SMINATOMA, SMINATOMB, SMINK, SMINKINC,
 31:      &                      SMINDISTSTART, SMINDISTFINISH, natomsina, natomsinb, natomsinc, atomsinalist, atomsinblist, 31:      &                      SMINDISTSTART, SMINDISTFINISH, natomsina, natomsinb, natomsinc, atomsinalist, atomsinblist,
 32:      &                      atomsinclist, atomsinalistlogical, atomsinblistlogical, atomsinclistlogical, ligcartstep, 32:      &                      atomsinclist, atomsinalistlogical, atomsinblistlogical, atomsinclistlogical, ligcartstep,
 33:      &                      ligtransstep, ligmovefreq, amchnmax, amchnmin, amchpmax, amchpmin, rotamert, rotmaxchange, 33:      &                      ligtransstep, ligmovefreq, amchnmax, amchnmin, amchpmax, amchpmin, rotamert, rotmaxchange,
 34:      &                      rotcentre, rotpselect, rotoccuw, rotcutoff, setchiralgeneric, PRMTOP, IGB, RGBMAX, CUT, 34:      &                      rotcentre, rotpselect, rotoccuw, rotcutoff, setchiralgeneric, PRMTOP, IGB, RGBMAX, CUT,
 35:      &                      SALTCON, macroiont, nmacroions, macroiondist 35:      &                      SALTCON, macroiont, nmacroions, macroiondist
 36:       USE modamber 36:       USE modamber
 37:       USE AMBER12_MUTATIONS, ONLY : AMBERMUTATION_SETUP!, AMBERMUTRB_SETUP 37:       USE AMBER12_MUTATIONS, ONLY : AMBERMUTATION_SETUP, AMBERMUTRB_SETUP
 38:       USE PORFUNCS 38:       USE PORFUNCS
 39:       USE MYGA_PARAMS 39:       USE MYGA_PARAMS
 40:       USE BGUPMOD 40:       USE BGUPMOD
 41:       USE GLJYMOD 41:       USE GLJYMOD
 42:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L 42:       USE CHIRO_MODULE, ONLY: CHIRO_SIGMA, CHIRO_MU, CHIRO_GAMMA, CHIRO_L
 43:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP 43:       USE CONVEX_POLYHEDRA_MODULE, ONLY: INITIALISE_POLYHEDRA, K_COMPRESS, K_OVERLAP
 44:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS, 44:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_MODE, LJ_GAUSS_RCUT, LJ_GAUSS_EPS,
 45:      &                        LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_PARAMS, 45:      &                        LJ_GAUSS_R0, LJ_GAUSS_SIGMASQ, LJ_GAUSS_PARAMS,
 46:      &                        LJ_GAUSS_INITIALISE 46:      &                        LJ_GAUSS_INITIALISE
 47:       USE OPP_MOD, ONLY: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS, 47:       USE OPP_MOD, ONLY: OPP_MODE, OPP_RCUT, OPP_K, OPP_PHI, OPP_PARAMS,
2157:         ELSE2157:         ELSE
2158:            WRITE(MYUNIT,'(A)') 'keyword> option not available for scoring mutation'2158:            WRITE(MYUNIT,'(A)') 'keyword> option not available for scoring mutation'
2159:            STOP2159:            STOP
2160:         ENDIF2160:         ENDIF
2161: 2161: 
2162:       ELSE IF (WORD.EQ.'AMBERMUTRIGID') THEN2162:       ELSE IF (WORD.EQ.'AMBERMUTRIGID') THEN
2163:         AMBERMUTRIGIDT = .TRUE.2163:         AMBERMUTRIGIDT = .TRUE.
2164:         RIGIDINIT = .TRUE.2164:         RIGIDINIT = .TRUE.
2165:         ATOMRIGIDCOORDT = .TRUE.2165:         ATOMRIGIDCOORDT = .TRUE.
2166:         AACONVERGENCET = .TRUE.2166:         AACONVERGENCET = .TRUE.
2167:         !CALL AMBERMUTRB_SETUP()2167:         CALL AMBERMUTRB_SETUP()
2168: 2168: 
2169:       ELSE IF (WORD.EQ.'AMBER9') THEN2169:       ELSE IF (WORD.EQ.'AMBER9') THEN
2170:         AMBERT=.TRUE.2170:         AMBERT=.TRUE.
2171:         WRITE(MYUNIT,'(A)') 'keyword> RADIUS set to 999 for AMBER9 run'2171:         WRITE(MYUNIT,'(A)') 'keyword> RADIUS set to 999 for AMBER9 run'
2172:         RADIUS=9992172:         RADIUS=999
2173: 2173: 
2174: !2174: !
2175: ! csw34> if residues are frozen with FREEZERES, call the amber routine2175: ! csw34> if residues are frozen with FREEZERES, call the amber routine
2176: ! to fill the FROZEN array correctly (in amberinterface.f)2176: ! to fill the FROZEN array correctly (in amberinterface.f)
2177: !2177: !


r33360/rbodyconfig_grouprotations.py 2017-09-29 15:30:10.190184713 +0100 r33359/rbodyconfig_grouprotations.py 2017-09-29 15:30:10.890193838 +0100
  1: import sys  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/AMBER/BHmutation_steps/rbodyconfig_grouprotations.py' in revision 33359
  2:  
  3: class protein(): 
  4:     def __init__(self, atom_dic, res_dic): 
  5:         self.atom_dic = atom_dic  # dictionary of all atoms with x,y,z coordinates and res id and name and atom name 
  6:         self.res_dic = res_dic  # dictionary containing a list of atom ids for all residues 
  7:  
  8:     def num_atoms(self):  # number of atoms 
  9:         return len(self.atom_dic) 
 10:  
 11:     def num_res(self):  # number of residues 
 12:         return len(self.res_dic) 
 13:  
 14:     # returns the name of a given residue 
 15:     def get_residue_name(self, residue): 
 16:         return atom_dic[res_dic[residue][0]][1] 
 17:  
 18:     # atom id from atom name given a residue id 
 19:     def get_atom_id_from_name(self, atom_name, res_num): 
 20:         for i in self.res_dic[res_num]: 
 21:             if self.atom_dic[i][0] == atom_name: 
 22:                 return i 
 23:  
 24:     # returns a list of atoms in a residue 
 25:     def get_atom_list(self, residue): 
 26:         all_atoms = self.atom_dic.keys() 
 27:         atomlist = [] 
 28:         for i in range(1, len(self.atom_dic) + 1): 
 29:             atom = self.atom_dic[i] 
 30:             if residue == atom[2]: 
 31:                 atomlist.append(all_atoms[i - 1]) 
 32:         return atomlist 
 33:  
 34: ################################################# 
 35: def reading_pdb(pdb): 
 36:     dic = {} 
 37:     atom_list = [] 
 38:     with open(pdb, 'r') as f_pdb: 
 39:         for line in f_pdb: 
 40:             if line.split()[0] == 'ATOM': 
 41:                 atom_list.append((line[12:16].strip(), line[17:20].strip(), int(line[22:26]))) 
 42:     for i in range(1, len(atom_list) + 1): 
 43:         dic[i] = atom_list[i - 1] 
 44:     return dic 
 45:  
 46: def reading_rbodyconfig(fname) 
 47:     dic = dict() 
 48:     with open(fname , "r") as f: 
 49:         for line in f: 
 50:             if line.split()[0] == 'GROUP': 
 51:                 dic[len(dic) + 1] = list() 
 52:             else: 
 53:                 dic[len(dic)].append(int(line.split()[0])) 
 54:     return dic 
 55:  
 56: def create_res_dic(atom_dic): 
 57:     res_dic = {} 
 58:     res_list = [] 
 59:     for i in atom_dic.keys(): 
 60:         res_list.append(atom_dic[i][2]) 
 61:     res_list = list(set(res_list)) 
 62:     for j in res_list: 
 63:         res_dic[j] = [] 
 64:     for k in atom_dic.keys(): 
 65:         res_dic[atom_dic[k][2]].append(k) 
 66:     return res_dic 
 67:  
 68: #################################################        
 69:  
 70: prob_uni = 0.025 
 71:  
 72: rbody_old = sys.argv[1] 
 73: shift = int(sys.argv[2]) 
 74: start = int(sys.argv[3]) 
 75: pdb_file = sys.argv[4] 
 76: atom_dic = reading_pdb(pdb_file) 
 77: res_dic = create_res_dic(atom_dic) 
 78: loaded_mol = protein(atom_dic, res_dic) 
 79:  
 80: amp_dic = {'ALA': (1.0,), 
 81:            'ASN': (0.5, 1.0), 
 82:            'ARG': (0.2, 0.3, 0.5), 
 83:            'ASP': (0.5, 1.0, 0.0), 
 84:            'CYS': (1.0,), 
 85:            'GLN': (0.3, 0.5, 1.0), 
 86:            'GLU': (0.3, 0.5, 1.0), 
 87:            'HIS': (0.3, 0.5), 
 88:            'HIE': (0.3, 0.5), 
 89:            'HIP': (0.3, 0.5), 
 90:            'HID': (0.3, 0.5), 
 91:            'ILE': (0.5, 1.0), 
 92:            'LEU': (0.5, 1.0), 
 93:            'LYS': (0.2, 0.3, 0.5), 
 94:            'MET': (0.5, 0.7), 
 95:            'PHE': (0.3, 0.5), 
 96:            'SER': (1.0,), 
 97:            'THR': (1.0,), 
 98:            'TRP': (0.3, 0.4), 
 99:            'TYR': (0.3, 0.5), 
100:            'VAL': (1.0,), 
101:            'AMB0': (0.1,)  # default for all non residue specific rotations 
102:           } 
103:  
104: #first read in old rbodyconfig and adjust it to give new output 
105: rb_dict = reading_rbodyconfig(rbody_old) 
106: used_atoms = list() 
107: for key in rb_dict.keys(): 
108:     if rb_dict[key][0] > start: 
109:         rb_dict[key] = [(atom + shift) for atom in rb_dict[key]] 
110:     used_atoms = used_atoms + rb_dict[key] 
111:  
112: rboutput = open("rbodyconfig.new" , "w") 
113: for key in rb_dict.keys(): 
114:     rboutput.write("GROUP " + str(len(rb_dict[key])) + "\n") 
115:     for atom in rb_dict[key]: 
116:         rboutput.write(str(atom) + "\n") 
117: rboutput.close() 
118:  
119: # number: [name, ax_atom1, ax_atom2, natom, amp, prob, [atom1, atom2, ...]] 
120: grot_dic = {} 
121: residues = res_dic.keys() 
122: for res in residues: 
123:     gr_res_name = loaded_mol.get_residue_name(res) 
124:     if gr_res_name not in ['PRO' , 'HYP']: 
125:         #N-CA rotation 
126:         prob = prob_uni 
127:         amp = 0.1  # use default amplitude 
128:  
129:         ax1 = loaded_mol.get_atom_id_from_name('N', res) 
130:         ax2 = loaded_mol.get_atom_id_from_name('CA', res) 
131:         gr_name = 'NCA' + str(res) 
132:         gr_atoms = [] 
133:         for gr_atom in loaded_mol.get_atom_list(res): 
134:             if atom_dic[gr_atom][0] not in ['N',  'C', 'CA', 'HA', 'O', 'OXT', 'H1', 'H2', 'H3']: 
135:                 gr_atoms.append(gr_atom) 
136:         inrbody = False 
137:         for atom in gr_atoms: 
138:             if atom in used_atoms: 
139:                 inrbody = True 
140:         if not(inrbody): 
141:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
142:  
143:     if gr_res_name not in ['PRO' , 'HYP']: 
144:         #C-CA rotation 
145:         prob = prob_uni 
146:         amp = 0.1  # use default amplitude 
147:         ax1 = loaded_mol.get_atom_id_from_name('C', res) 
148:         ax2 = loaded_mol.get_atom_id_from_name('CA', res) 
149:         gr_name = 'CCA' + str(res) 
150:         gr_atoms = [] 
151:         for gr_atom in loaded_mol.get_atom_list(res): 
152:             if atom_dic[gr_atom][0] not in ['N',  'C', 'CA', 'HA', 'O', 'OXT', 'H1', 'H2', 'H3']: 
153:                 gr_atoms.append(gr_atom) 
154:         inrbody = False 
155:         for atom in gr_atoms: 
156:             if atom in used_atoms: 
157:                 inrbody = True 
158:         if not(inrbody): 
159:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
160:  
161:     if gr_res_name not in ['PRO', 'HYP' , 'ALA', 'GLY']: 
162:         # CA-CB rotation 
163:         prob = prob_uni 
164:         try: 
165:             amp = amp_dic[gr_res_name][0]  # use default amplitude 
166:         except KeyError: 
167:             amp = 0.1 
168:         ax1 = loaded_mol.get_atom_id_from_name('CA', res) 
169:         ax2 = loaded_mol.get_atom_id_from_name('CB', res) 
170:         gr_name = 'CACB' + str(res) 
171:         gr_atoms = [] 
172:         for gr_atom in loaded_mol.get_atom_list(res): 
173:             if atom_dic[gr_atom][0] not in ['N', 'H', 'C', 'CA', 'HA', 'CB', 'O', 'OXT', 'H1', 'H2', 'H3']: 
174:                 gr_atoms.append(gr_atom) 
175:         inrbody = False 
176:         for atom in gr_atoms: 
177:             if atom in used_atoms: 
178:                 inrbody = True 
179:         if not(inrbody): 
180:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
181:  
182:     if gr_res_name not in ['PRO', 'HYP' , 'ILE', 'ALA', 'GLY', 'SER', 'CYS', 'THR', 'VAL']: 
183:         # CB-CG rotation 
184:         prob = prob_uni 
185:         try: 
186:             amp = amp_dic[gr_res_name][1]  # use default amplitude 
187:         except KeyError: 
188:             amp = 0.1 
189:         ax1 = loaded_mol.get_atom_id_from_name('CB', res) 
190:         ax2 = loaded_mol.get_atom_id_from_name('CG', res) 
191:         gr_name = 'CBCG' + str(res) 
192:         gr_atoms = [] 
193:         for gr_atom in loaded_mol.get_atom_list(res): 
194:             if atom_dic[gr_atom][0] not in ['N', 'H', 'C', 'CA', 'HA', 'CB', 'CG', 'HB', 'HB1', 'HB2', 
195:                                             'HB3', 'O', 'OXT', 'H1', 'H2', 'H3']: 
196:                 gr_atoms.append(gr_atom) 
197:         inrbody = False 
198:         for atom in gr_atoms: 
199:             if atom in used_atoms: 
200:                 inrbody = True 
201:         if not(inrbody): 
202:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
203:  
204:     if gr_res_name in ['ARG', 'LYS', 'GLU', 'GLN']: 
205:         # CG-CD rotation 
206:         prob = prob_uni 
207:         try: 
208:             amp = amp_dic[gr_res_name][2]  # use default amplitude 
209:         except KeyError: 
210:             amp = 0.1 
211:         ax1 = loaded_mol.get_atom_id_from_name('CG', res) 
212:         ax2 = loaded_mol.get_atom_id_from_name('CD', res) 
213:         gr_name = 'CGCD' + str(res) 
214:         gr_atoms = [] 
215:         atom_exc = ['N', 'H', 'C', 'CA', 'HA', 'CB', 'CG', 'HB', 'HB1', 'HB2', 'HB3', 'O', 'OXT', 'H1', 
216:                     'H2', 'H3', 'CD', 'HG1', 'HG2', 'HG11', 'HG12', 'HG13', 'HG21', 'HG22', 'HG23'] 
217:         for gr_atom in loaded_mol.get_atom_list(res): 
218:             if atom_dic[gr_atom][0] not in atom_exc: 
219:                 gr_atoms.append(gr_atom) 
220:                 inrbody = False 
221:         for atom in gr_atoms: 
222:             if atom in used_atoms: 
223:                 inrbody = True 
224:         if not(inrbody): 
225:             grot_dic[len(grot_dic) + 1] = [gr_name, ax1, ax2, len(gr_atoms), amp, prob, gr_atoms] 
226:  
227:    
228: # write output for group rotation files 
229: groups_file = open('atomgroups', 'w') 
230: for group in range(1, len(grot_dic.keys()) + 1): 
231:     string = 'GROUP ' + grot_dic[group][0] + ' ' + str(grot_dic[group][1]) + ' ' + str(grot_dic[group][2]) 
232:     string += ' ' + str(grot_dic[group][3]) + ' ' + str(grot_dic[group][4]) + ' ' + str(grot_dic[group][5]) + "\n" 
233:     for atom in grot_dic[group][6]: 
234:         string += str(atom) + "\n" 
235:     groups_file.write(string) 
236: groups_file.close() 
237:  
238:  
239:  


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0