hdiff output

r33459/amber_mutations.F90 2017-11-10 13:30:17.522462614 +0000 r33458/amber_mutations.F90 2017-11-10 13:30:17.774465976 +0000
 26:   DOUBLE PRECISION, ALLOCATABLE, SAVE :: CA_REFERENCE(:) 26:   DOUBLE PRECISION, ALLOCATABLE, SAVE :: CA_REFERENCE(:)
 27:  27: 
 28:   CONTAINS 28:   CONTAINS
 29:   !setup mutational system, initialise coordinates correctly, and according to the right size 29:   !setup mutational system, initialise coordinates correctly, and according to the right size
 30:   SUBROUTINE AMBERMUTATION_SETUP()   30:   SUBROUTINE AMBERMUTATION_SETUP()  
 31:      IMPLICIT NONE 31:      IMPLICIT NONE
 32:      INTEGER :: MIUNIT,GETUNIT,J1,J2,NENTRIES,NTERMINI,TESTINT 32:      INTEGER :: MIUNIT,GETUNIT,J1,J2,NENTRIES,NTERMINI,TESTINT
 33:      LOGICAL :: YESNO , NTERT 33:      LOGICAL :: YESNO , NTERT
 34:      CHARACTER(200) :: ENTRY_ 34:      CHARACTER(200) :: ENTRY_
 35:      CHARACTER(25) , DIMENSION(:) , ALLOCATABLE :: ENTRIES 35:      CHARACTER(25) , DIMENSION(:) , ALLOCATABLE :: ENTRIES
 36:      INTEGER, PARAMETER :: NMUTOCC = 1000  
 37:  36: 
 38:      !check there is a file contianing the mutational information 37:      !check there is a file contianing the mutational information
 39:      YESNO = .FALSE. 38:      YESNO = .FALSE.
 40:      INQUIRE(FILE='amber_mutations',EXIST=YESNO) 39:      INQUIRE(FILE='amber_mutations',EXIST=YESNO)
 41:      IF (.NOT.YESNO) THEN 40:      IF (.NOT.YESNO) THEN
 42:         WRITE(MYUNIT,'(A)') ' ambermut> No mutation information given' 41:         WRITE(MYUNIT,'(A)') ' ambermut> No mutation information given'
 43:         STOP 42:         STOP
 44:      ENDIF 43:      ENDIF
 45:      !get the number of residues, and their atom positions 44:      !get the number of residues, and their atom positions
 46:      CALL TOPOLOGY_READER() 45:      CALL TOPOLOGY_READER()
 47:      !open the mutation information 46:      !open the mutation information
 48:      MIUNIT = GETUNIT() 47:      MIUNIT = GETUNIT()
 49:      OPEN(UNIT=MIUNIT,FILE='amber_mutations',status='unknown') 48:      OPEN(UNIT=MIUNIT,FILE='amber_mutations',status='unknown')
 50:      WRITE(MYUNIT,*) 'ambermut> Reading in mutations allowed' 49:      WRITE(MYUNIT,*) 'ambermut> Reading in mutations allowed'
 51:      READ(MIUNIT,*) NRESMUT 50:      READ(MIUNIT,*) NRESMUT
 52:      WRITE(MYUNIT,*) 'ambermut> ',NRESIDUES,' residues, of which ',NRESMUT,'can be mutated' 51:      WRITE(MYUNIT,*) 'ambermut> ',NRESIDUES,' residues, of which ',NRESMUT,'can be mutated'
 53:      !all allocations and reallocations are taken care of here, as we will not chnage the number of residues 52:      !all allocations and reallocations are taken care of here, as we will not chnage the number of residues
  53:      ALLOCATE(TERMINI_RES(NRESIDUES))
 54:      IF (ALLOCATED(MUTATION_INFO)) DEALLOCATE(MUTATION_INFO) 54:      IF (ALLOCATED(MUTATION_INFO)) DEALLOCATE(MUTATION_INFO)
 55:      ALLOCATE(MUTATION_INFO(NRESMUT)) 55:      ALLOCATE(MUTATION_INFO(NRESMUT))
 56:      IF (ALLOCATED(PREVIOUS_MUTATION)) DEALLOCATE(PREVIOUS_MUTATION) 56:      IF (ALLOCATED(PREVIOUS_MUTATION)) DEALLOCATE(PREVIOUS_MUTATION)
 57:      ALLOCATE(PREVIOUS_MUTATION(NRESMUT)) 57:      ALLOCATE(PREVIOUS_MUTATION(NRESMUT))
 58:      IF (ALLOCATED(TERMINI_RES)) DEALLOCATE(TERMINI_RES) 58:      IF (ALLOCATED(TERMINI_RES)) DEALLOCATE(TERMINI_RES)
 59:      ALLOCATE(TERMINI_RES(NRESIDUES)) 59:      ALLOCATE(TERMINI_RES(NRESIDUES))
 60:      TERMINI_RES(:) = 0 60:      TERMINI_RES(:) = 0
 61:      !best sequence 
 62:      IF (ALLOCATED(BESTMUTSEQ)) DEALLOCATE(BESTMUTSEQ) 
 63:      ALLOCATE(BESTMUTSEQ(NRESIDUES)) 
 64:      BESTMUTSEQ(:) = "    " 
 65:      !all sequences including rejected ones 
 66:      IF (ALLOCATED(MUTSEQ)) DEALLOCATE(MUTSEQ) 
 67:      ALLOCATE(MUTSEQ(NMUTOCC,NRESIDUES)) 
 68:      MUTSEQ(:,:) = "    " 
 69:      !storing information about rejecting sequences 
 70:      IF (ALLOCATED(SEQREJECTED)) DEALLOCATE(SEQREJECTED) 
 71:      ALLOCATE(SEQREJECTED(NMUTOCC)) 
 72:      SEQREJECTED(:) = .FALSE. 
 73:      NSEQSTORED = 0 
 74:      !first sequence to be stored is initial configuration 
 75:      MUTSEQ(1,:) = AMBER12_RESNAME(:) 
 76:      NSEQSTORED = NSEQSTORED + 1 
 77:      !read next line, this contains the terminal residues 61:      !read next line, this contains the terminal residues
 78:      READ(MIUNIT,*) NTERMINI 62:      READ(MIUNIT,*) NTERMINI
 79:      READ(MIUNIT,'(A)',END=101) ENTRY_               !read line 63:      READ(MIUNIT,'(A)',END=101) ENTRY_               !read line
 80:      ALLOCATE(ENTRIES(NTERMINI)) 64:      ALLOCATE(ENTRIES(NTERMINI))
 81:      ENTRIES(:)='' 65:      ENTRIES(:)=''
 82:      CALL READ_LINE(ENTRY_,NTERMINI,ENTRIES) 66:      CALL READ_LINE(ENTRY_,NTERMINI,ENTRIES)
 83:      !entries now contains the information about the termini, we need to differentiate between C and N termini 67:      !entries now contains the information about the termini, we need to differentiate between C and N termini
 84:      !the first entry ought to be the N terminus and then we switch between then, this might cause problems for ACE, NME, NHE!!! 68:      !the first entry ought to be the N terminus and then we switch between then, this might cause problems for ACE, NME, NHE!!!
 85:      NTERT = .TRUE. 69:      NTERT = .TRUE.
 86:      DO J1=1,NRESIDUES 70:      DO J1=1,NRESIDUES
178:      RETURN162:      RETURN
179:   END SUBROUTINE AMBERMUTRB_SETUP163:   END SUBROUTINE AMBERMUTRB_SETUP
180:   164:   
181:   !mutate protein165:   !mutate protein
182:   SUBROUTINE AMBERMUT_STEP(COORDINATES , RESNUMBER)166:   SUBROUTINE AMBERMUT_STEP(COORDINATES , RESNUMBER)
183:      INTEGER , INTENT(OUT) :: RESNUMBER167:      INTEGER , INTENT(OUT) :: RESNUMBER
184:      CHARACTER(LEN=4) :: OLDRES , OLDRES1 , NEWRES , NEWRES1168:      CHARACTER(LEN=4) :: OLDRES , OLDRES1 , NEWRES , NEWRES1
185:      CHARACTER(LEN=6) :: NMUT_STRING , STARTINDEX_STRING, START_STR, SHIFT_STR169:      CHARACTER(LEN=6) :: NMUT_STRING , STARTINDEX_STRING, START_STR, SHIFT_STR
186:      CHARACTER(LEN=50) :: OPTION_STRING, PATH_STRING170:      CHARACTER(LEN=50) :: OPTION_STRING, PATH_STRING
187:      DOUBLE PRECISION :: COORDINATES(3*NATOMS)171:      DOUBLE PRECISION :: COORDINATES(3*NATOMS)
188:      INTEGER ::  SHIFT, START, J1172:      INTEGER ::  SHIFT, START
189:  173:  
190:      !let's store all information first in case we have to go back!174:      !let's store all information first in case we have to go back!
191:      PREVIOUS_MUTATION = MUTATION_INFO175:      PREVIOUS_MUTATION = MUTATION_INFO
192:      !we have a new mutation176:      !we have a new mutation
193:      NMUTATION = NMUTATION + 1177:      NMUTATION = NMUTATION + 1
194:      WRITE(NMUT_STRING,'(I6)') NMUTATION - 1 178:      WRITE(NMUT_STRING,'(I6)') NMUTATION - 1 
195:      !before we do anything, we save the old lowest minima179:      !before we do anything, we save the old lowest minima
196:      CALL AMBERMUT_CURR_LOWEST()180:      CALL AMBERMUT_CURR_LOWEST()
197:      !select a residue to mutate181:      !select a residue to mutate
198:      CALL SELECT_MUTATION(RESNUMBER , OLDRES1 , NEWRES1)182:      CALL SELECT_MUTATION(RESNUMBER , OLDRES1 , NEWRES1)
199:      !if it is a terminal residue, we need to go for a different set of atoms and coordinates in the coordinate creation script183:      !if it is a terminal residue, we need to go for a different set of atoms and coordinates in the coordinate creation script
200:      IF (TERMINI_RES(RESNUMBER).EQ.1) THEN184:      IF (TERMINI_RES(RESNUMBER).EQ.1) THEN
201:         OLDRES = "C" // OLDRES1185:         OLDRES = "C" // OLDRES1
202:         NEWRES = "C" // NEWRES1186:         NEWRES = "C" // NEWRES1
203:      ELSE IF (TERMINI_RES(RESNUMBER).EQ.2) THEN187:      ELSE IF (TERMINI_RES(RESNUMBER).EQ.2) THEN
204:         OLDRES = "N" // OLDRES1188:         OLDRES = "N" // OLDRES1
205:         NEWRES = "N" // NEWRES1189:         NEWRES = "N" // NEWRES1
206:      ELSE190:      ELSE
207:         OLDRES = OLDRES1191:         OLDRES = OLDRES1
208:         NEWRES = NEWRES1192:         NEWRES = NEWRES1
209:      ENDIF193:      ENDIF                
210:      WRITE(MUTUNIT,'(A)') 'Currently the best match to the objective function is:' 
211:      WRITE(MUTUNIT,'(A,F20.10)') 'Mutation score: ' , BESTMUTSCORE 
212:      WRITE(MUTUNIT,'(A)',ADVANCE='NO') 'Sequence: ' 
213:      DO J1=1,NRESIDUES-1 
214:         WRITE(MUTUNIT,'(A)',ADVANCE='NO') BESTMUTSEQ(J1) // "  " 
215:      ENDDO 
216:      WRITE(MUTUNIT,'(A)') BESTMUTSEQ(NRESIDUES) 
217:      WRITE(MUTUNIT,'(A)') '==============================' 
218:      WRITE(MUTUNIT,'(A,I6,4A)') 'Mutate residue ' , RESNUMBER , ' from ' , OLDRES , ' to ' , NEWRES194:      WRITE(MUTUNIT,'(A,I6,4A)') 'Mutate residue ' , RESNUMBER , ' from ' , OLDRES , ' to ' , NEWRES
219:      WRITE(STARTINDEX_STRING,'(I6)') AMBER12_RESSTART(RESNUMBER)195:      WRITE(STARTINDEX_STRING,'(I6)') AMBER12_RESSTART(RESNUMBER)
220: 196: 
221:      !dump the coordinates for the old residue, and move things to safety197:      !dump the coordinates for the old residue, and move things to safety
222:      CALL DUMP_RESIDUE_COORDS(RESNUMBER , COORDINATES)198:      CALL DUMP_RESIDUE_COORDS(RESNUMBER , COORDINATES)
223:      CALL SYSTEM('mv coords.prmtop coords.prmtop.'//TRIM(ADJUSTL(NMUT_STRING)))199:      CALL SYSTEM('mv coords.prmtop coords.prmtop.'//TRIM(ADJUSTL(NMUT_STRING)))
224:      CALL SYSTEM('mv coords.inpcrd coords.inpcrd.'//TRIM(ADJUSTL(NMUT_STRING)))200:      CALL SYSTEM('mv coords.inpcrd coords.inpcrd.'//TRIM(ADJUSTL(NMUT_STRING)))
225:      CALL SYSTEM('mv start start.'//TRIM(ADJUSTL(NMUT_STRING)))201:      CALL SYSTEM('mv start start.'//TRIM(ADJUSTL(NMUT_STRING)))
226:      CALL SYSTEM('mv atomgroups atomgroups.'//TRIM(ADJUSTL(NMUT_STRING)))202:      CALL SYSTEM('mv atomgroups atomgroups.'//TRIM(ADJUSTL(NMUT_STRING)))
227:      !create mutated coordinates and a new perm.allow file203:      !create mutated coordinates and a new perm.allow file
258: #endif234: #endif
259:         CALL SYSTEM('mv rbodyconfig.new rbodyconfig.')235:         CALL SYSTEM('mv rbodyconfig.new rbodyconfig.')
260:         CALL SYSTEM('mv coordsinirigid coordsinirigid.'//TRIM(ADJUSTL(NMUT_STRING)))236:         CALL SYSTEM('mv coordsinirigid coordsinirigid.'//TRIM(ADJUSTL(NMUT_STRING)))
261:         CALL SYSTEM('cp start coordsinirigid')237:         CALL SYSTEM('cp start coordsinirigid')
262:      ENDIF238:      ENDIF
263:      CALL SYSTEM('rm tmp.pdb')239:      CALL SYSTEM('rm tmp.pdb')
264:      !finally reinitialise AMBER with new groups, coordinates and topology240:      !finally reinitialise AMBER with new groups, coordinates and topology
265:      CALL REINITIALISE_AMBER()241:      CALL REINITIALISE_AMBER()
266:      !now remove old chiral states used for checking (the rest is done when we initialise the chirality in mc.F)242:      !now remove old chiral states used for checking (the rest is done when we initialise the chirality in mc.F)
267:      CALL DEALLOC_STATES_MUTATION()243:      CALL DEALLOC_STATES_MUTATION()
268:      !finally store the new sequence244:      RETURN
269:      NSEQSTORED = NSEQSTORED + 1 
270:      MUTSEQ(NSEQSTORED, :) = AMBER12_RESNAME(:) 
271:      RETURN      
272:   END SUBROUTINE AMBERMUT_STEP245:   END SUBROUTINE AMBERMUT_STEP
273: 246: 
274:   SUBROUTINE SELECT_MUTATION(RESNUMBER , OLDRES , NEWRES)247:   SUBROUTINE SELECT_MUTATION(RESNUMBER , OLDRES , NEWRES)
275:      INTEGER , INTENT(OUT) :: RESNUMBER248:      INTEGER , INTENT(OUT) :: RESNUMBER
276:      CHARACTER(LEN=4) , INTENT(OUT) :: OLDRES , NEWRES 249:      CHARACTER(LEN=4) , INTENT(OUT) :: OLDRES , NEWRES 
277:      CHARACTER(LEN=4) :: SELECTED_MUT250:      CHARACTER(LEN=4) :: SELECTED_MUT
278:      INTEGER :: ENTRIES , NCURR , J1 , SELECTED_ID , SELECTED_RES251:      INTEGER :: ENTRIES , NCURR , J1 , SELECTED_ID , SELECTED_RES
279:      DOUBLE PRECISION :: PROB_RES_SELECT(NRESMUT,2) , NMUTATED , PROB , PROBTOT , RANDOM, DPRAND252:      DOUBLE PRECISION :: PROB_RES_SELECT(NRESMUT,2) , NMUTATED , PROB , PROBTOT , RANDOM, DPRAND
280:      DOUBLE PRECISION , ALLOCATABLE :: PROB_MUT_SELECT(:,:)253:      DOUBLE PRECISION , ALLOCATABLE :: PROB_MUT_SELECT(:,:)
281:      !create probability array to select residue id254:      !create probability array to select residue id
599:      CALL SYSTEM('cp atomgroups.'//TRIM(ADJUSTL(NMUT_STRING))//' atomgroups')572:      CALL SYSTEM('cp atomgroups.'//TRIM(ADJUSTL(NMUT_STRING))//' atomgroups')
600:      CALL SYSTEM('cp perm.allow.'//TRIM(ADJUSTL(NMUT_STRING))//' perm.allow')573:      CALL SYSTEM('cp perm.allow.'//TRIM(ADJUSTL(NMUT_STRING))//' perm.allow')
601:      IF (AMBERMUTRIGIDT) THEN574:      IF (AMBERMUTRIGIDT) THEN
602:         CALL SYSTEM('cp rbodyconfig.'//TRIM(ADJUSTL(NMUT_STRING))//' rbodyconfig') 575:         CALL SYSTEM('cp rbodyconfig.'//TRIM(ADJUSTL(NMUT_STRING))//' rbodyconfig') 
603:         CALL SYSTEM('cp start coordsinirigid')     576:         CALL SYSTEM('cp start coordsinirigid')     
604:      ENDIF577:      ENDIF
605:      !now reinitialise once more578:      !now reinitialise once more
606:      CALL REINITIALISE_AMBER()579:      CALL REINITIALISE_AMBER()
607:      !reset chirality580:      !reset chirality
608:      CALL DEALLOC_STATES_MUTATION()581:      CALL DEALLOC_STATES_MUTATION()
609:      !store new (old) sequence and note the previous one was rejected 
610:      SEQREJECTED(NSEQSTORED) = .TRUE. 
611:      NSEQSTORED = NSEQSTORED + 1 
612:      MUTSEQ(NSEQSTORED, :) = AMBER12_RESNAME(:) 
613:      RETURN 
614:   END SUBROUTINE REVERSE_MUTATION582:   END SUBROUTINE REVERSE_MUTATION
615: 583: 
616:   !scoring function to judge how good mutation is584:   !scoring function to judge how good mutation is
617:   SUBROUTINE MUTATION_E(SCORE,COORDS,MODE,TERMID)585:   SUBROUTINE MUTATION_E(SCORE,COORDS,MODE,TERMID)
618:      DOUBLE PRECISION, INTENT(OUT) :: SCORE586:      DOUBLE PRECISION, INTENT(OUT) :: SCORE
619:      DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)587:      DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)
620:      INTEGER, INTENT(IN) :: MODE, TERMID588:      INTEGER, INTENT(IN) :: MODE, TERMID
621:      DOUBLE PRECISION :: DPRAND, EREAL, GRADATOMS(3*NATOMS)589:      DOUBLE PRECISION :: DPRAND, EREAL, GRADATOMS(3*NATOMS)
622:      INTEGER :: ATOMID, PARMEDUNIT, GETUNIT, J1590:      INTEGER :: ATOMID, PARMEDUNIT, GETUNIT, J1
623:      TYPE(POT_ENE_REC_C) :: DECOMPOSED_E591:      TYPE(POT_ENE_REC_C) :: DECOMPOSED_E


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0