hdiff output

r33119/amber_mutations.F90 2017-08-04 18:30:25.069023139 +0100 r33118/amber_mutations.F90 2017-08-04 18:30:25.961035005 +0100
119:         ENDDO119:         ENDDO
120:      ENDDO120:      ENDDO
121: 101  CONTINUE121: 101  CONTINUE
122:      CLOSE(MIUNIT)122:      CLOSE(MIUNIT)
123:      !call the grouprotation set up here (not in keywords)123:      !call the grouprotation set up here (not in keywords)
124:      CALL MUT_SETUP_GROUPROTATION(1,.FALSE.,.FALSE.,0)         124:      CALL MUT_SETUP_GROUPROTATION(1,.FALSE.,.FALSE.,0)         
125:      RETURN125:      RETURN
126:   END SUBROUTINE AMBERMUTATION_SETUP126:   END SUBROUTINE AMBERMUTATION_SETUP
127:   127:   
128:   !mutate protein128:   !mutate protein
129:   SUBROUTINE AMBERMUT_STEP(COORDINATES , RESNUMBER)129:   SUBROUTINE AMBERMUT_STEP(COORDINATES)
130:      INTEGER , INTENT(OUT) :: RESNUMBER130:      INTEGER :: RESNUMBER
131:      CHARACTER(LEN=4) :: OLDRES , OLDRES1 , NEWRES , NEWRES1131:      CHARACTER(LEN=4) :: OLDRES , OLDRES1 , NEWRES , NEWRES1
132:      CHARACTER(LEN=6) :: NMUT_STRING , STARTINDEX_STRING132:      CHARACTER(LEN=6) :: NMUT_STRING , STARTINDEX_STRING
133:      CHARACTER(LEN=25) :: OPTION_STRING133:      CHARACTER(LEN=25) :: OPTION_STRING
134:      DOUBLE PRECISION :: COORDINATES(3*NATOMS)134:      DOUBLE PRECISION :: COORDINATES(3*NATOMS)
135:  135:  
136:      !let's store all information first in case we have to go back!136:      !let's store all information first in case we have to go back!
137:      PREVIOUS_MUTATION = MUTATION_INFO137:      PREVIOUS_MUTATION = MUTATION_INFO
138:      !we have a new mutation138:      !we have a new mutation
139:      NMUTATION = NMUTATION + 1139:      NMUTATION = NMUTATION + 1
140:      WRITE(NMUT_STRING,'(I6)') NMUTATION - 1 140:      WRITE(NMUT_STRING,'(I6)') NMUTATION - 1 
457:         QMINNATOMS(1:NSAVE)=NATOMS ! to prevent reading from uninitialised memory457:         QMINNATOMS(1:NSAVE)=NATOMS ! to prevent reading from uninitialised memory
458:         COORDSO(1:3*NATOMS,1:NPAR)=0.0D0458:         COORDSO(1:3*NATOMS,1:NPAR)=0.0D0
459:         VT(1:NATOMS)=0.0D0459:         VT(1:NATOMS)=0.0D0
460:         VAT(1:NATOMS,1:NPAR)=0.0D0460:         VAT(1:NATOMS,1:NPAR)=0.0D0
461:         DO J1=1,NSAVE461:         DO J1=1,NSAVE
462:            QMIN(J1)=1.0D10462:            QMIN(J1)=1.0D10
463:            NPCALL_QMIN(J1)=0463:            NPCALL_QMIN(J1)=0
464:         ENDDO464:         ENDDO
465:   END SUBROUTINE REINITIALISE_AMBER465:   END SUBROUTINE REINITIALISE_AMBER
466: 466: 
467:   SUBROUTINE REVERSE_MUTATION(RESNUMBER)467:   SUBROUTINE REVERSE_MUTATION()
468:      CHARACTER(LEN=6) :: NMUT_STRING468:      CHARACTER(LEN=6) :: NMUT_STRING     
469:      INTEGER :: STARTATOM, FINALATOM_OLD,FINALATOM_NEW,SHIFT 
470:      INTEGER , INTENT(IN) :: RESNUMBER     
471:      469:      
472:      !reload the correct information into MUTATION_INFO470:      !reload the correct information into MUTATION_INFO
473:      MUTATION_INFO = PREVIOUS_MUTATION 471:      MUTATION_INFO = PREVIOUS_MUTATION 
474:      WRITE(NMUT_STRING,'(I6)') NMUTATION - 1472:      WRITE(NMUT_STRING,'(I6)') NMUTATION - 1
475:      STARTATOM = AMBER12_RESSTART(RESNUMBER) 
476:      FINALATOM_OLD = AMBER12_RESEND(RESNUMBER) 
477:      !move all the files we need back into place (we use the lowest previous minimum to restart)473:      !move all the files we need back into place (we use the lowest previous minimum to restart)
478:      CALL SYSTEM('cp coords.prmtop.'//TRIM(ADJUSTL(NMUT_STRING))//' coords.prmtop')474:      CALL SYSTEM('cp coords.prmtop.'//TRIM(ADJUSTL(NMUT_STRING))//' coords.prmtop')
479:      CALL TOPOLOGY_READER() 
480:      !get the changed number of atoms 
481:      FINALATOM_NEW = AMBER12_RESEND(RESNUMBER) 
482:      SHIFT = FINALATOM_NEW - FINALATOM_OLD 
483:      !correct wrong information (we havent reinitialised yet, so NATOMS is still wrong) 
484:      AMBER12_RESEND(NRESIDUES) = AMBER12_RESEND(NRESIDUES) + SHIFT 
485:      AMBER12_RESNATOM(NRESIDUES) = AMBER12_RESNATOM(NRESIDUES) + SHIFT 
486:      !create final input files needed 
487:      CALL SYSTEM('cp coords.'//TRIM(ADJUSTL(NMUT_STRING))//'.1.rst coords.inpcrd')475:      CALL SYSTEM('cp coords.'//TRIM(ADJUSTL(NMUT_STRING))//'.1.rst coords.inpcrd')
488:      CALL SYSTEM('cp start.'//TRIM(ADJUSTL(NMUT_STRING))//' start')  !this one is the wrong file?476:      CALL SYSTEM('cp start.'//TRIM(ADJUSTL(NMUT_STRING))//' start')  !this one is the wrong file?
489:      CALL SYSTEM('cp atomgroups.'//TRIM(ADJUSTL(NMUT_STRING))//' atomgroups')477:      CALL SYSTEM('cp atomgroups.'//TRIM(ADJUSTL(NMUT_STRING))//' ambergroups')
490:      CALL SYSTEM('cp perm.allow.'//TRIM(ADJUSTL(NMUT_STRING))//' perm.allow')478:      CALL SYSTEM('cp perm.allow.'//TRIM(ADJUSTL(NMUT_STRING))//' perm.allow')
491:      !now reinitialise once more479:      !now reinitialise once more
492:      CALL REINITIALISE_AMBER()480:      CALL REINITIALISE_AMBER()
493:   END SUBROUTINE REVERSE_MUTATION481:   END SUBROUTINE REVERSE_MUTATION
494: 482: 
495:   !scoring function to judge how good mutation is483:   !scoring function to judge how good mutation is
496:   SUBROUTINE MUTATION_E(SCORE,COORDS,MODE,TERMID)484:   SUBROUTINE MUTATION_E(SCORE,COORDS)
497:      DOUBLE PRECISION, INTENT(OUT) :: SCORE485:      DOUBLE PRECISION, INTENT(OUT) :: SCORE
498:      DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)486:      DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)     
499:      INTEGER, INTENT(IN) :: MODE, TERMID487:      DOUBLE PRECISION :: DPRAND
500:      DOUBLE PRECISION :: DPRAND, EREAL, GRADATOMS(3*NATOMS) 
501:      INTEGER :: ATOMID, PARMEDUNIT, GETUNIT, J1 
502:      TYPE(POT_ENE_REC_C) :: DECOMPOSED_E 
503:      CHARACTER(200) ENTRY_ 
504:      INTEGER , PARAMETER :: NWORDS=20 
505:      CHARACTER(25) :: ENTRIES(NWORDS)='' 
506:      CHARACTER(LEN=6) :: J1_STRING 
507:    488:    
508:      IF (MODE.EQ.1) THEN489:      SCORE=DPRAND()
509:         SCORE=DPRAND() 
510:      ELSE IF (MODE.EQ.2) THEN 
511:         CALL AMBER12_ENERGY_AND_GRADIENT(NATOMS, COORDS, EREAL, GRADATOMS, DECOMPOSED_E) 
512:         WRITE(MUTUNIT,'(A)') 'Energy decomposition' 
513:         WRITE(MUTUNIT,'(A,F20.10)') 'Total energy:        ', DECOMPOSED_E % TOTAL 
514:         WRITE(MUTUNIT,'(A,F20.10)') 'Total van der Waals: ', DECOMPOSED_E % VDW_TOT 
515:         WRITE(MUTUNIT,'(A,F20.10)') 'Total electronic:    ', DECOMPOSED_E % ELEC_TOT 
516:         WRITE(MUTUNIT,'(A,F20.10)') 'Generalised Born:    ', DECOMPOSED_E % GB 
517:         WRITE(MUTUNIT,'(A,F20.10)') 'Surface energy:      ', DECOMPOSED_E % SURF 
518:         WRITE(MUTUNIT,'(A,F20.10)') 'Bond energy:         ', DECOMPOSED_E % BOND 
519:         WRITE(MUTUNIT,'(A,F20.10)') 'Angular term:        ', DECOMPOSED_E % ANGLE 
520:         WRITE(MUTUNIT,'(A,F20.10)') 'Dihedral term:       ', DECOMPOSED_E % DIHEDRAL 
521:         WRITE(MUTUNIT,'(A,F20.10)') 'vdW 1-4 term:        ', DECOMPOSED_E % VDW_14 
522:         WRITE(MUTUNIT,'(A,F20.10)') 'Electronic 1-4:      ', DECOMPOSED_E % ELEC_14 
523:         WRITE(MUTUNIT,'(A,F20.10)') 'Restraints:          ', DECOMPOSED_E % RESTRAINT 
524:         WRITE(MUTUNIT,'(A,F20.10)') 'Urey Bradley angle:  ', DECOMPOSED_E % ANGLE_UB 
525:         WRITE(MUTUNIT,'(A,F20.10)') 'Improper energy:     ', DECOMPOSED_E % IMP 
526:         WRITE(MUTUNIT,'(A,F20.10)') 'CMAP:                ', DECOMPOSED_E % CMAP 
527:         IF (TERMID.EQ.0) THEN 
528:            SCORE = DECOMPOSED_E % TOTAL 
529:         ELSE IF (TERMID.EQ.1) THEN 
530:            SCORE = DECOMPOSED_E % VDW_TOT 
531:         ELSE IF (TERMID.EQ.2) THEN 
532:            SCORE = DECOMPOSED_E % ELEC_TOT 
533:         ELSE IF (TERMID.EQ.3) THEN 
534:            SCORE = DECOMPOSED_E % GB 
535:         ELSE IF (TERMID.EQ.4) THEN 
536:            SCORE = DECOMPOSED_E % SURF 
537:         ELSE IF (TERMID.EQ.5) THEN 
538:            SCORE = DECOMPOSED_E % BOND 
539:         ELSE IF (TERMID.EQ.6) THEN 
540:            SCORE = DECOMPOSED_E % ANGLE 
541:         ELSE IF (TERMID.EQ.7) THEN 
542:            SCORE = DECOMPOSED_E % DIHEDRAL 
543:         ELSE IF (TERMID.EQ.8) THEN 
544:            SCORE = DECOMPOSED_E % VDW_14 
545:         ELSE IF (TERMID.EQ.9) THEN 
546:            SCORE = DECOMPOSED_E % ELEC_14 
547:         ELSE IF (TERMID.EQ.10) THEN 
548:            SCORE = DECOMPOSED_E % RESTRAINT 
549:         ELSE IF (TERMID.EQ.11) THEN 
550:            SCORE = DECOMPOSED_E % ANGLE_UB 
551:         ELSE IF (TERMID.EQ.12) THEN 
552:            SCORE = DECOMPOSED_E % IMP 
553:         ELSE IF (TERMID.EQ.13) THEN 
554:            SCORE = DECOMPOSED_E % CMAP 
555:         ENDIF 
556:      ELSE IF (MODE.EQ.3) THEN 
557:         ATOMID = AMBER12_RESEND(TERMID) !get last atom in first group 
558:         !open new file to write parmed input 
559:         PARMEDUNIT = GETUNIT() 
560:         OPEN(PARMEDUNIT,FILE='parmed_in',STATUS='NEW') 
561:         WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') 'addExclusions @1'         
562:         DO J1=2,ATOMID 
563:            WRITE(J1_STRING,'(I6)') J1 
564:            WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING)) 
565:         ENDDO 
566:         WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ' @1'        
567:         DO J1=2,ATOMID-1 
568:            WRITE(J1_STRING,'(I6)') J1 
569:            WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING)) 
570:         ENDDO 
571:         WRITE(J1_STRING,'(I6)') ATOMID 
572:         WRITE(PARMEDUNIT,'(A)') ','//TRIM(ADJUSTL(J1_STRING)) 
573:         WRITE(J1_STRING,'(I6)') ATOMID+1 
574:         WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') 'addExclusions @'//TRIM(ADJUSTL(J1_STRING))         
575:         DO J1=ATOMID+2,NATOMS 
576:            WRITE(J1_STRING,'(I6)') J1 
577:            WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING)) 
578:         ENDDO 
579:         WRITE(J1_STRING,'(I6)') ATOMID+1 
580:         WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ' @'//TRIM(ADJUSTL(J1_STRING))        
581:         DO J1=ATOMID+2,NATOMS-1 
582:            WRITE(J1_STRING,'(I6)') J1 
583:            WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING)) 
584:         ENDDO 
585:         WRITE(J1_STRING,'(I6)') NATOMS-1 
586:         WRITE(PARMEDUNIT,'(A)') ','//TRIM(ADJUSTL(J1_STRING)) 
587:         WRITE(PARMEDUNIT,'(A)') 'loadRestrt current.inpcrd' 
588:         WRITE(J1_STRING,'(I6)') AMBERMUTIGB 
589:         WRITE(PARMEDUNIT,'(A)') 'energy cutoff 15.0 igb '//TRIM(ADJUSTL(J1_STRING))//' saltcon 0.1' 
590:         WRITE(PARMEDUNIT,'(A)') 'quit' 
591:         CLOSE(PARMEDUNIT) 
592:         CALL AMBER12_WRITE_RESTART(COORDS, 'current.inpcrd',LEN('current.inpcrd')) 
593:         !create new topology without interactions and calculate energy 
594:         CALL SYSTEM('parmed.py -n coords.prmtop parmed_in > parmed_out') 
595:         OPEN(PARMEDUNIT,FILE='parmed_out',STATUS='OLD') 
596:         DO 
597:           ENTRIES(:)='' 
598:           READ(PARMEDUNIT,'(A)',END=588) ENTRY_ 
599:           CALL READ_LINE(ENTRY_,NWORDS,ENTRIES) 
600:           IF (ENTRIES(1).EQ.'TOTAL') THEN 
601:              READ(ENTRIES(3),'(F20.10)') SCORE 
602:              GOTO 588 
603:           ENDIF 
604:         ENDDO 
605: 588     CONTINUE 
606:         CALL SYSTEM('rm current.inpcrd parmed_in parmed_out')       
607:      ENDIF 
608:   END SUBROUTINE MUTATION_E490:   END SUBROUTINE MUTATION_E
609: 491: 
610:   SUBROUTINE AMBERMUT_CURR_LOWEST()492:   SUBROUTINE AMBERMUT_CURR_LOWEST()
611:      CHARACTER(LEN=6) :: J1_STRING, NMUT_STRING493:      CHARACTER(LEN=6) :: J1_STRING, NMUT_STRING
612:      INTEGER :: J1    494:      INTEGER :: J1    
613:  495:  
614:      !to save the previous best strutures 496:      !to save the previous best strutures 
615:      DO J1=1,NSAVE497:      DO J1=1,NSAVE
616:         WRITE(J1_STRING,'(I6)') J1498:         WRITE(J1_STRING,'(I6)') J1
617:         WRITE(NMUT_STRING,'(I6)') NMUTATION - 1499:         WRITE(NMUT_STRING,'(I6)') NMUTATION - 1


r33119/commons.f90 2017-08-04 18:30:25.289026065 +0100 r33118/commons.f90 2017-08-04 18:30:26.181037929 +0100
654:       INTEGER, ALLOCATABLE ::  MLQOUTCOME(:)654:       INTEGER, ALLOCATABLE ::  MLQOUTCOME(:)
655:       INTEGER, ALLOCATABLE ::  LJADDNN(:,:)655:       INTEGER, ALLOCATABLE ::  LJADDNN(:,:)
656: 656: 
657:       INTEGER, DIMENSION(:,:), ALLOCATABLE :: BONDS !for QCIAMBER657:       INTEGER, DIMENSION(:,:), ALLOCATABLE :: BONDS !for QCIAMBER
658: 658: 
659: !OPEP interface659: !OPEP interface
660:       LOGICAL :: OPEPT, OPEP_RNAT660:       LOGICAL :: OPEPT, OPEP_RNAT
661: 661: 
662: !AMBER mutational steps662: !AMBER mutational steps
663:       LOGICAL :: AMBERMUTATIONT663:       LOGICAL :: AMBERMUTATIONT
664:       INTEGER :: MUTUNIT,NMUTATION,MUTATIONFREQ,MUTTESTSTEPS,AMBERMUTFF,AMBERMUTIGB,MUTENERGY,MUTTERMID664:       INTEGER :: MUTUNIT,NMUTATION,MUTATIONFREQ,MUTTESTSTEPS,AMBERMUTFF,AMBERMUTIGB
665: 665: 
666: !Orbital variables666: !Orbital variables
667:       LOGICAL :: ORBITALS667:       LOGICAL :: ORBITALS
668:       INTEGER :: NROTS, NORBS, ORBVAREXPONENT668:       INTEGER :: NROTS, NORBS, ORBVAREXPONENT
669:       DOUBLE PRECISION, ALLOCATABLE :: R2INTS(:,:), DIPINTS(:,:,:)669:       DOUBLE PRECISION, ALLOCATABLE :: R2INTS(:,:), DIPINTS(:,:,:)
670: END MODULE COMMONS670: END MODULE COMMONS


r33119/keywords.f 2017-08-04 18:30:25.521029150 +0100 r33118/keywords.f 2017-08-04 18:30:26.409040960 +0100
1238:       ORBVAREXPONENT = -1!default to easily identifiable (quasi-)nonsense value.1238:       ORBVAREXPONENT = -1!default to easily identifiable (quasi-)nonsense value.
1239: 1239: 
1240: ! AMBER mutations1240: ! AMBER mutations
1241:       AMBERMUTATIONT = .FALSE.1241:       AMBERMUTATIONT = .FALSE.
1242:       MUTATIONFREQ = 10001242:       MUTATIONFREQ = 1000
1243:       MUTTESTSTEPS = 501243:       MUTTESTSTEPS = 50
1244:       NMUTATION = 01244:       NMUTATION = 0
1245:       MUTUNIT = GETUNIT()1245:       MUTUNIT = GETUNIT()
1246:       AMBERMUTIGB = 21246:       AMBERMUTIGB = 2
1247:       AMBERMUTFF = 141247:       AMBERMUTFF = 14
1248:       MUTENERGY = 1 
1249:       MUTTERMID = 0 
1250: 1248: 
1251:       CALL FILE_OPEN('data', DATA_UNIT, .FALSE.)1249:       CALL FILE_OPEN('data', DATA_UNIT, .FALSE.)
1252: 1250: 
1253: !      OPEN (5,FILE='data',STATUS='OLD')1251: !      OPEN (5,FILE='data',STATUS='OLD')
1254: 1252: 
1255: !190   CALL INPUT(END,5)1253: !190   CALL INPUT(END,5)
1256: 190   CALL INPUT(END, DATA_UNIT)1254: 190   CALL INPUT(END, DATA_UNIT)
1257:       IF (.NOT. END) THEN1255:       IF (.NOT. END) THEN
1258:         CALL READU(WORD)1256:         CALL READU(WORD)
1259:       ENDIF1257:       ENDIF
2097:       ELSE IF (WORD.EQ.'AMBERMUTFF') THEN2095:       ELSE IF (WORD.EQ.'AMBERMUTFF') THEN
2098:         CALL READI(AMBERMUTFF)2096:         CALL READI(AMBERMUTFF)
2099:         IF (AMBERMUTFF.EQ.14) THEN2097:         IF (AMBERMUTFF.EQ.14) THEN
2100:            WRITE(MYUNIT,'(A)') 'keyword> use ff14SB force field for mutations'2098:            WRITE(MYUNIT,'(A)') 'keyword> use ff14SB force field for mutations'
2101:         ELSE IF (AMBERMUTFF.EQ.99) THEN2099:         ELSE IF (AMBERMUTFF.EQ.99) THEN
2102:            WRITE(MYUNIT,'(A)') 'keyword> use ff99SB force field for mutations'2100:            WRITE(MYUNIT,'(A)') 'keyword> use ff99SB force field for mutations'
2103:         ELSE2101:         ELSE
2104:            WRITE(MYUNIT,'(A)') 'keyword> invalid force field choice'2102:            WRITE(MYUNIT,'(A)') 'keyword> invalid force field choice'
2105:            STOP2103:            STOP
2106:         ENDIF2104:         ENDIF
2107:       2105:        
2108:       ELSE IF (WORD.EQ.'AMBERMUTENERGY') THEN2106: 
2109:         CALL READI(MUTENERGY) 
2110:         CALL READI(MUTTERMID) 
2111:         IF (MUTENERGY.EQ.1) THEN 
2112:            WRITE(MYUNIT,'(A)') 'keyword> use random numbers to score mutations' 
2113:         ELSE IF (MUTENERGY.EQ.2) THEN 
2114:            WRITE(MYUNIT,'(A)') 'keyword> use decomposed energy to score mutations' 
2115:         ELSE IF (MUTENERGY.EQ.3) THEN 
2116:            WRITE(MYUNIT,'(A)') 'keyword> use interaction energy between parts to score mutations' 
2117:            WRITE(MYUNIT,'(A,I6)') 'keyword> first component ends at residue ' , MUTTERMID 
2118:         ELSE 
2119:            WRITE(MYUNIT,'(A)') 'keyword> option not available for scoring mutation' 
2120:            STOP 
2121:         ENDIF 
2122: 2107: 
2123:       ELSE IF (WORD.EQ.'AMBER9') THEN2108:       ELSE IF (WORD.EQ.'AMBER9') THEN
2124:         AMBERT=.TRUE.2109:         AMBERT=.TRUE.
2125:         WRITE(MYUNIT,'(A)') 'keyword> RADIUS set to 999 for AMBER9 run'2110:         WRITE(MYUNIT,'(A)') 'keyword> RADIUS set to 999 for AMBER9 run'
2126:         RADIUS=9992111:         RADIUS=999
2127: 2112: 
2128: !2113: !
2129: ! csw34> if residues are frozen with FREEZERES, call the amber routine2114: ! csw34> if residues are frozen with FREEZERES, call the amber routine
2130: ! to fill the FROZEN array correctly (in amberinterface.f)2115: ! to fill the FROZEN array correctly (in amberinterface.f)
2131: !2116: !


r33119/mc.F 2017-08-04 18:30:25.741032077 +0100 r33118/mc.F 2017-08-04 18:30:26.629043887 +0100
 97:  97: 
 98: ! ds656 > energy before and after BGUPTA swap (without relaxation) 98: ! ds656 > energy before and after BGUPTA swap (without relaxation)
 99:       DOUBLE PRECISION :: EBSWAP=0.0D0, EASWAP=0.0D0, EASWAPQ=0.0D0, DE1=0.0D0, DE2=0.0D0 99:       DOUBLE PRECISION :: EBSWAP=0.0D0, EASWAP=0.0D0, EASWAPQ=0.0D0, DE1=0.0D0, DE2=0.0D0
100: 100: 
101:       DOUBLE PRECISION QTEMP(3,NATOMSALLOC), Q4, Q6101:       DOUBLE PRECISION QTEMP(3,NATOMSALLOC), Q4, Q6
102: 102: 
103: ! kr366> AMBER mutation steps103: ! kr366> AMBER mutation steps
104:       DOUBLE PRECISION :: EMUTP = 0.0D0 , EMUTN=0.0D0104:       DOUBLE PRECISION :: EMUTP = 0.0D0 , EMUTN=0.0D0
105:       LOGICAL :: MUTATEDT=.FALSE. , DOMUTATIONSTEPT=.FALSE.105:       LOGICAL :: MUTATEDT=.FALSE. , DOMUTATIONSTEPT=.FALSE.
106:       CHARACTER(LEN=10) :: NMUT_STRING106:       CHARACTER(LEN=10) :: NMUT_STRING
107:       INTEGER :: MUTATEDRES = 0 
108: 107: 
109: !kr366> as we may need to change the size of the coordinate arrays for mutational steps they are now allocatable108: !kr366> as we may need to change the size of the coordinate arrays for mutational steps they are now allocatable
110:       ALLOCATE(BESTCOORDS(3*NATOMSALLOC,NPAR))109:       ALLOCATE(BESTCOORDS(3*NATOMSALLOC,NPAR))
111:       ALLOCATE(SAVECOORDS(3*NATOMSALLOC))110:       ALLOCATE(SAVECOORDS(3*NATOMSALLOC))
112:       ALLOCATE(TEMPCOORDS(3*NATOMSALLOC))111:       ALLOCATE(TEMPCOORDS(3*NATOMSALLOC))
113:       ALLOCATE(SCREENC(3*NATOMSALLOC))112:       ALLOCATE(SCREENC(3*NATOMSALLOC))
114: 113: 
115:       QNEWRES=0114:       QNEWRES=0
116:       QGCBH=0115:       QGCBH=0
117: 116: 
1200:                IF (MACROIONT) THEN1199:                IF (MACROIONT) THEN
1201:                 CALL MACROION_MOVES(J1,COORDS(:,JP),movableatomlist,nmovableatoms,ligmovet,blockmovet,nblocks,1200:                 CALL MACROION_MOVES(J1,COORDS(:,JP),movableatomlist,nmovableatoms,ligmovet,blockmovet,nblocks,
1202:      1                              atomsinblock,STEP(JP))1201:      1                              atomsinblock,STEP(JP))
1203:                END IF1202:                END IF
1204: ! normal AMBER12 step taking (without MACROION model)              1203: ! normal AMBER12 step taking (without MACROION model)              
1205:                ELSE IF (AMBER12T.AND..NOT.MACROIONT) THEN1204:                ELSE IF (AMBER12T.AND..NOT.MACROIONT) THEN
1206:                   !kr366> mutational step,  otherwise continue with the normal steps 1205:                   !kr366> mutational step,  otherwise continue with the normal steps 
1207:                   IF (AMBERMUTATIONT.AND.DOMUTATIONSTEPT) THEN1206:                   IF (AMBERMUTATIONT.AND.DOMUTATIONSTEPT) THEN
1208:                      WRITE(MUTUNIT,'(A,I8)') 'Number of atoms before mutation: ',NATOMS1207:                      WRITE(MUTUNIT,'(A,I8)') 'Number of atoms before mutation: ',NATOMS
1209:                      WRITE(MYUNIT,'(A,I8)') ' mc> Number of atoms before mutation: ',NATOMS1208:                      WRITE(MYUNIT,'(A,I8)') ' mc> Number of atoms before mutation: ',NATOMS
1210:                      CALL MUTATION_E(EMUTP,COORDS(:,JP),MUTENERGY,MUTTERMID) !get the energy for the old molecule before we reinitialise AMBER!1209:                      CALL MUTATION_E(EMUTP,COORDS(:,JP)) !get the energy for the old molecule before we reinitialise AMBER!
1211:                      CALL AMBERMUT_STEP(COORDS(:,JP),MUTATEDRES)1210:                      CALL AMBERMUT_STEP(COORDS(:,JP))
1212:                      IF (SETCHIRAL) THEN1211:                      IF (SETCHIRAL) THEN
1213:                         WRITE(MYUNIT,'(A)') ' mc> Storing chiral information for initial (mutated) structure'1212:                         WRITE(MYUNIT,'(A)') ' mc> Storing chiral information for initial (mutated) structure'
1214:                         CALL INIT_CHIRAL(COORDS(:,1))1213:                         CALL INIT_CHIRAL(COORDS(:,1))
1215:                      END IF1214:                      END IF
1216:                      IF (NOCISTRANS) THEN1215:                      IF (NOCISTRANS) THEN
1217:                         WRITE(MYUNIT,'(A)') ' mc> Storing cis/trans information for initial (mutated) structure'1216:                         WRITE(MYUNIT,'(A)') ' mc> Storing cis/trans information for initial (mutated) structure'
1218:                         CALL INIT_CIS_TRANS(COORDS(:,1))1217:                         CALL INIT_CIS_TRANS(COORDS(:,1))
1219:                      ENDIF1218:                      ENDIF
1220:                      WRITE(MUTUNIT,'(A,I8)') 'Number of atoms after mutation: ',NATOMS  1219:                      WRITE(MUTUNIT,'(A,I8)') 'Number of atoms after mutation: ',NATOMS  
1221:                      WRITE(MYUNIT,'(A,I8)') ' mc> Number of atoms after mutation: ',NATOMS                    1220:                      WRITE(MYUNIT,'(A,I8)') ' mc> Number of atoms after mutation: ',NATOMS                    
2308:                   POTEL=GCEBEST(JP)2307:                   POTEL=GCEBEST(JP)
2309:                   EPREV(JP)=GCEBEST(JP)2308:                   EPREV(JP)=GCEBEST(JP)
2310:                   IF (DEBUG) WRITE(MYUNIT,'(A,G20.10)') 'mc> potel and eprev reset to ',GCEBEST(JP)2309:                   IF (DEBUG) WRITE(MYUNIT,'(A,G20.10)') 'mc> potel and eprev reset to ',GCEBEST(JP)
2311:                   IF (DEBUG) WRITE(MYUNIT,'(A,I6)') 'mc> number of atoms reset to ',NATOMS2310:                   IF (DEBUG) WRITE(MYUNIT,'(A,I6)') 'mc> number of atoms reset to ',NATOMS
2312:                ENDIF2311:                ENDIF
2313:             ENDIF2312:             ENDIF
2314: 2313: 
2315:             !the mutational steps are accepted for better bonding properties, but the energies used 2314:             !the mutational steps are accepted for better bonding properties, but the energies used 
2316:             !otherwise refer to the best potential energies for one sequence2315:             !otherwise refer to the best potential energies for one sequence
2317:             IF (MUTATEDT.AND.(MOD((J1-MUTTESTSTEPS),MUTATIONFREQ).EQ.0)) THEN !kr366> here we check our mutational step               2316:             IF (MUTATEDT.AND.(MOD((J1-MUTTESTSTEPS),MUTATIONFREQ).EQ.0)) THEN !kr366> here we check our mutational step               
2318:                CALL MUTATION_E(EMUTN,TEMPCOORDS,MUTENERGY,MUTTERMID) !get the energy for the new molecule2317:                CALL MUTATION_E(EMUTN,TEMPCOORDS) !get the energy for the new molecule
2319:                WRITE(MYUNIT, '(2(A,F15.8))') ' mc> Energy before mutation: ',EMUTP,' | energy after mutation: ',EMUTN2318:                WRITE(MYUNIT, '(2(A,F15.8))') ' mc> Energy before mutation: ',EMUTP,' | energy after mutation: ',EMUTN
2320:                WRITE(MUTUNIT, '(2(A,F15.8))') 'Energy before mutation: ',EMUTP,' | energy after mutation: ',EMUTN2319:                WRITE(MUTUNIT, '(2(A,F15.8))') 'Energy before mutation: ',EMUTP,' | energy after mutation: ',EMUTN
2321:                IF (EMUTN.LT.EMUTP) THEN2320:                IF (EMUTN.LT.EMUTP) THEN
2322:                   RANDOM=0.0D02321:                   RANDOM=0.0D0
2323:                   ATEST=.TRUE.2322:                   ATEST=.TRUE.
2324:                ELSE2323:                ELSE
2325:                   RANDOM=DPRAND()2324:                   RANDOM=DPRAND()
2326:                   DUMMY = EMUTN - EMUTP2325:                   DUMMY = EMUTN - EMUTP
2327:                   DUMMY=EXP(-DUMMY/MAX(MCTEMP,1.0D-100))2326:                   DUMMY=EXP(-DUMMY/MAX(MCTEMP,1.0D-100))
2328:                   IF (DUMMY.GT.RANDOM) THEN2327:                   IF (DUMMY.GT.RANDOM) THEN
2329:                      ATEST=.TRUE.2328:                      ATEST=.TRUE.
2330:                   ELSE2329:                   ELSE
2331:                      ATEST=.FALSE.2330:                      ATEST=.FALSE.
2332:                   ENDIF2331:                   ENDIF
2333:                ENDIF2332:                ENDIF
2334:                IF (ATEST) THEN2333:                IF (ATEST) THEN
2335:                   WRITE(MUTUNIT,'(A)') 'Mutation accpeted'2334:                   WRITE(MUTUNIT,'(A)') 'Mutation accpeted'
2336:                   WRITE(MYUNIT, '(A)') ' mc> Mutation accepted'2335:                   WRITE(MYUNIT, '(A)') ' mc> Mutation accepted'
2337:                ELSE2336:                ELSE
2338:                   CALL REVERSE_MUTATION(MUTATEDRES)2337:                   CALL REVERSE_MUTATION()
2339:                   IF (SETCHIRAL) THEN2338:                   IF (SETCHIRAL) THEN
2340:                      WRITE(MYUNIT,'(A)') ' mc> Storing chiral information for initial (mutated) structure'2339:                      WRITE(MYUNIT,'(A)') ' mc> Storing chiral information for initial (mutated) structure'
2341:                      CALL INIT_CHIRAL(COORDS(:,1))2340:                      CALL INIT_CHIRAL(COORDS(:,1))
2342:                   END IF2341:                   END IF
2343:                   IF (NOCISTRANS) THEN2342:                   IF (NOCISTRANS) THEN
2344:                      WRITE(MYUNIT,'(A)') ' mc> Storing cis/trans information for initial (mutated) structure'2343:                      WRITE(MYUNIT,'(A)') ' mc> Storing cis/trans information for initial (mutated) structure'
2345:                      CALL INIT_CIS_TRANS(COORDS(:,1))2344:                      CALL INIT_CIS_TRANS(COORDS(:,1))
2346:                   ENDIF2345:                   ENDIF
2347:                   WRITE(MYUNIT, '(A)')  ' mc> Calculating initial energy after reversing mutation'2346:                   WRITE(MYUNIT, '(A)')  ' mc> Calculating initial energy after reversing mutation'
2348:                   EPSSAVE=EPSSPHERE2347:                   EPSSAVE=EPSSPHERE


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0