hdiff output

r33512/amber_mutations.F90 2017-11-22 11:30:15.369987743 +0000 r33511/amber_mutations.F90 2017-11-22 11:30:15.957995555 +0000
1344:       ENDDO1344:       ENDDO
1345:   END SUBROUTINE READ_LINE1345:   END SUBROUTINE READ_LINE
1346:  1346:  
1347:   SUBROUTINE DISTANCE(A1,A2,DIST)1347:   SUBROUTINE DISTANCE(A1,A2,DIST)
1348:      DOUBLE PRECISION, INTENT(IN) :: A1(3), A2(3)1348:      DOUBLE PRECISION, INTENT(IN) :: A1(3), A2(3)
1349:      DOUBLE PRECISION, INTENT(OUT) :: DIST1349:      DOUBLE PRECISION, INTENT(OUT) :: DIST
1350: 1350: 
1351:      DIST = SQRT((A1(1)-A2(1))**2+(A1(2)-A2(2))**2+(A1(3)-A2(3))**2)1351:      DIST = SQRT((A1(1)-A2(1))**2+(A1(2)-A2(2))**2+(A1(3)-A2(3))**2)
1352:   END SUBROUTINE DISTANCE1352:   END SUBROUTINE DISTANCE
1353: 1353: 
1354:   SUBROUTINE AMBER12DECOMP(NA,X,DECOMPOSED) 
1355:      INTEGER, INTENT(IN) :: NA 
1356:      DOUBLE PRECISION, INTENT(IN) :: X(3*NA) 
1357:      DOUBLE PRECISION :: DUMMYE 
1358:      DOUBLE PRECISION :: DUMMYGRAD(3*NA) 
1359:      TYPE(POT_ENE_REC_C) ,INTENT(OUT) :: DECOMPOSED 
1360:  
1361:      CALL AMBER12_ENERGY_AND_GRADIENT(NATOMS, COORDS, DUMMYE, DUMMYGRAD, DECOMPOSED) 
1362:   END SUBROUTINE AMBER12DECOMP 
1363:  
1364: END MODULE AMBER12_MUTATIONS1354: END MODULE AMBER12_MUTATIONS


r33512/mc.F 2017-11-22 11:30:15.673991775 +0000 r33511/mc.F 2017-11-22 11:30:16.249999438 +0000
 26:  26: 
 27:       USE QMODULE , ONLY : QMIN, QMINP, INTEQMIN 27:       USE QMODULE , ONLY : QMIN, QMINP, INTEQMIN
 28:       USE modcharmm 28:       USE modcharmm
 29:       USE MODAMBER9, ONLY : MDSTEPT,CISARRAY1,CHIARRAY1,NOCISTRANSDNA,NOCISTRANSRNA, 29:       USE MODAMBER9, ONLY : MDSTEPT,CISARRAY1,CHIARRAY1,NOCISTRANSDNA,NOCISTRANSRNA,
 30:      &                      SETCHIRAL,AMCHPMAX,DOLIGMOVE,LIGMOVEFREQ, AMCHNMAX, ROTAMERT, 30:      &                      SETCHIRAL,AMCHPMAX,DOLIGMOVE,LIGMOVEFREQ, AMCHNMAX, ROTAMERT,
 31:      &                      E_IGB, E_BOND, E_ANGLE, E_DIHEDRAL, 31:      &                      E_IGB, E_BOND, E_ANGLE, E_DIHEDRAL,
 32:      &                      E_VDW, E_14_VDW, E_ELEC, E_14_ELEC, IGB, MACROIONT 32:      &                      E_VDW, E_14_VDW, E_ELEC, E_14_ELEC, IGB, MACROIONT
 33:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, POT_ENE_REC_C 33:       USE AMBER12_INTERFACE_MOD, ONLY : AMBER12_WRITE_RESTART, AMBER12_WRITE_PDB, POT_ENE_REC_C
 34:       USE CHIRALITY, ONLY : INIT_CHIRAL, INIT_CIS_TRANS 34:       USE CHIRALITY, ONLY : INIT_CHIRAL, INIT_CIS_TRANS
 35:       USE AMBER12_MUTATIONS, ONLY: AMBERMUT_STEP, REVERSE_MUTATION , MUTATION_E, PRINT_CURRENT_SEQ, 35:       USE AMBER12_MUTATIONS, ONLY: AMBERMUT_STEP, REVERSE_MUTATION , MUTATION_E, PRINT_CURRENT_SEQ,
 36:      &                             AMBERMUTDUMP, AMBER12DECOMP 36:      &                             AMBERMUTDUMP
 37:       USE porfuncs 37:       USE porfuncs
 38:       USE AMHGLOBALS, ONLY: NMRES,OMOVI,AVEP,NUMPRO,IRES 38:       USE AMHGLOBALS, ONLY: NMRES,OMOVI,AVEP,NUMPRO,IRES
 39:       USE AMH_INTERFACES, ONLY:E_WRITE 39:       USE AMH_INTERFACES, ONLY:E_WRITE
 40:       USE ROTAMER 40:       USE ROTAMER
 41:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_TAKESTEP 41:       USE LJ_GAUSS_MOD, ONLY: LJ_GAUSS_TAKESTEP
 42:       USE OPP_MOD, ONLY: OPP_TAKESTEP 42:       USE OPP_MOD, ONLY: OPP_TAKESTEP
 43:       USE BOX_DERIVATIVES, ONLY: BD_TAKESTEP 43:       USE BOX_DERIVATIVES, ONLY: BD_TAKESTEP
 44:  44: 
 45:       IMPLICIT NONE 45:       IMPLICIT NONE
 46: #ifdef MPI 46: #ifdef MPI
 87:       COMMON /TOT/ NQTOT 87:       COMMON /TOT/ NQTOT
 88:       COMMON /Q4C/ QSTART, QFINISH 88:       COMMON /Q4C/ QSTART, QFINISH
 89:  89: 
 90:       character(len=10)       :: datechar,timechar,zonechar 90:       character(len=10)       :: datechar,timechar,zonechar
 91:       integer                 :: values(8),itime1 91:       integer                 :: values(8),itime1
 92:       DOUBLE PRECISION :: DISTGROUPX2,DISTGROUPY2,DISTGROUPZ2,DISTGROUPCENTRE,DUMMY 92:       DOUBLE PRECISION :: DISTGROUPX2,DISTGROUPY2,DISTGROUPZ2,DISTGROUPCENTRE,DUMMY
 93:       integer :: J6, L, M, SUMSQUAREDIFF, SUMSQUAREDIFF1D, HBONDDIFF 93:       integer :: J6, L, M, SUMSQUAREDIFF, SUMSQUAREDIFF1D, HBONDDIFF
 94:       INTEGER GETUNIT 94:       INTEGER GETUNIT
 95:       INTEGER DUMPUNIQUEUNIT 95:       INTEGER DUMPUNIQUEUNIT
 96:  96: 
 97: ! khs26> Energy decomposition for AMBER 12, getting it through the routines written in the mutational step interface (kr366) 97: ! khs26> Energy decomposition for AMBER 12
 98:       TYPE(POT_ENE_REC_C) :: A12DECOMPOSED_E 98: !      TYPE(POT_ENE_REC_C) :: AMBER12_ENERGY_DECOMP
 99:  99: 
100: ! ds656 > energy before and after BGUPTA swap (without relaxation)100: ! ds656 > energy before and after BGUPTA swap (without relaxation)
101:       DOUBLE PRECISION :: EBSWAP=0.0D0, EASWAP=0.0D0, EASWAPQ=0.0D0, DE1=0.0D0, DE2=0.0D0101:       DOUBLE PRECISION :: EBSWAP=0.0D0, EASWAP=0.0D0, EASWAPQ=0.0D0, DE1=0.0D0, DE2=0.0D0
102: 102: 
103:       DOUBLE PRECISION QTEMP(3,NATOMSALLOC), Q4, Q6103:       DOUBLE PRECISION QTEMP(3,NATOMSALLOC), Q4, Q6
104: 104: 
105: ! kr366> AMBER mutation steps105: ! kr366> AMBER mutation steps
106:       DOUBLE PRECISION :: EMUTP = 0.0D0 , EMUTN=0.0D0 , MUTSCORE=0.0D0 106:       DOUBLE PRECISION :: EMUTP = 0.0D0 , EMUTN=0.0D0 , MUTSCORE=0.0D0 
107:       LOGICAL :: MUTATEDT=.FALSE. , DOMUTATIONSTEPT=.FALSE.107:       LOGICAL :: MUTATEDT=.FALSE. , DOMUTATIONSTEPT=.FALSE.
108:       CHARACTER(LEN=10) :: NMUT_STRING108:       CHARACTER(LEN=10) :: NMUT_STRING
860: ! mo361> rigid body rotation move switch860: ! mo361> rigid body rotation move switch
861:                IF(TRANSLATERIGIDT.AND.MOD(J1,TRANSLATERIGIDFREQ).EQ.0) THEN861:                IF(TRANSLATERIGIDT.AND.MOD(J1,TRANSLATERIGIDFREQ).EQ.0) THEN
862:                        DOTRANSLATERIGID=.TRUE.862:                        DOTRANSLATERIGID=.TRUE.
863:                ENDIF863:                ENDIF
864: 864: 
865:                IF(ROTATEHINGET .AND. MOD(J1,ROTATEHINGEFREQ).EQ.0) THEN865:                IF(ROTATEHINGET .AND. MOD(J1,ROTATEHINGEFREQ).EQ.0) THEN
866:                        DOROTATEHINGE=.TRUE.866:                        DOROTATEHINGE=.TRUE.
867:                ENDIF867:                ENDIF
868: ! kr366> mutation moves for AMBER12               868: ! kr366> mutation moves for AMBER12               
869:                IF(AMBERMUTATIONT.AND.(MOD(J1,MUTATIONFREQ).EQ.0).AND.(J1.GT.0)) THEN869:                IF(AMBERMUTATIONT.AND.(MOD(J1,MUTATIONFREQ).EQ.0).AND.(J1.GT.0)) THEN
870:                        IF((NSTEPS-J1).GT.MUTTESTSTEPS) THEN !check we can fit the mutation plus exploration into the run870:                        WRITE(MUTUNIT,'(A,I10,A)') 'Mutational step after ',J1,' steps'
871:                           WRITE(MUTUNIT,'(A,I10,A)') 'Mutational step after ',J1,' steps'871:                        DOMUTATIONSTEPT=.TRUE.
872:                           DOMUTATIONSTEPT=.TRUE.872:                        MUTATEDT=.TRUE.        !as we need to consider different tests if we mutate!
873:                           MUTATEDT=.TRUE.        !as we need to consider different tests if we mutate!873:                        NMUTSTEP=J1
874:                           NMUTSTEP=J1 
875:                           DOGROUPROT=.FALSE. 
876:                        ELSE 
877:                           DOGROUPROT=.TRUE. 
878:                        ENDIF 
879:                ELSE IF (AMBERMUTATIONT.AND..NOT.MOD(J1,MUTATIONFREQ).EQ.0) THEN874:                ELSE IF (AMBERMUTATIONT.AND..NOT.MOD(J1,MUTATIONFREQ).EQ.0) THEN
880:                        DOGROUPROT=.TRUE.875:                        DOGROUPROT=.FALSE.
881:                ENDIF876:                ENDIF
882: 877: 
883: !878: !
884: ! csw34> Coordinates are saved so that moves can be undone879: ! csw34> Coordinates are saved so that moves can be undone
885: !880: !
886:                SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)881:                SAVECOORDS(1:3*NATOMS)=COORDS(1:3*NATOMS,JP)
887:                if (boxderivt) savebox_params(1:6) = box_params(1:6) ! dj337: save box params882:                if (boxderivt) savebox_params(1:6) = box_params(1:6) ! dj337: save box params
888: ! csw34> If you want to look at the effect of moves, you can dump out883: ! csw34> If you want to look at the effect of moves, you can dump out
889: ! the structure BEFORE the move here.884: ! the structure BEFORE the move here.
890: !                 CALL A9DUMPPDB(COORDS(:,JP),"beforemove")885: !                 CALL A9DUMPPDB(COORDS(:,JP),"beforemove")
1649:                         IF (IGB.EQ.0) THEN1644:                         IF (IGB.EQ.0) THEN
1650:                            WRITE(MYUNIT,'(A10,F20.10)') 'hbond:    ', E_IGB1645:                            WRITE(MYUNIT,'(A10,F20.10)') 'hbond:    ', E_IGB
1651: !     khs26> igb = 10 is the Poisson-Boltzmann case1646: !     khs26> igb = 10 is the Poisson-Boltzmann case
1652:                         ELSE IF (IGB.EQ.10) THEN1647:                         ELSE IF (IGB.EQ.10) THEN
1653:                            WRITE(MYUNIT,'(A10,F20.10)') 'ipb:      ', E_IGB1648:                            WRITE(MYUNIT,'(A10,F20.10)') 'ipb:      ', E_IGB
1654: !     khs26> other values of igb are detailed in the AMBER manual1649: !     khs26> other values of igb are detailed in the AMBER manual
1655:                         ELSE1650:                         ELSE
1656:                            WRITE(MYUNIT,'(A10,F20.10)') 'igb:      ', E_IGB1651:                            WRITE(MYUNIT,'(A10,F20.10)') 'igb:      ', E_IGB
1657:                         END IF1652:                         END IF
1658:                      ELSE IF (AMBER12T) THEN1653:                      ELSE IF (AMBER12T) THEN
1659:                         CALL AMBER12DECOMP(NATOMS,COORDS(:,JP),A12DECOMPOSED_E)1654: ! khs26> Add the energy decomposition terms from AMBER12 here! I need to think of how to get access to it from
1660:                         WRITE(MUTUNIT,'(A)') 'Energy decomposition'1655: ! potential.f without adding it to commons (eugh).
1661:                         WRITE(MUTUNIT,'(A,F20.10)') 'Total energy:        ', A12DECOMPOSED_E % TOTAL1656:                         WRITE(MYUNIT,'(A80)') 'Energy decomposition not yet implemented for AMB12GMIN.' 
1662:                         WRITE(MUTUNIT,'(A,F20.10)') 'Total van der Waals: ', A12DECOMPOSED_E % VDW_TOT 
1663:                         WRITE(MUTUNIT,'(A,F20.10)') 'Total electronic:    ', A12DECOMPOSED_E % ELEC_TOT 
1664:                         WRITE(MUTUNIT,'(A,F20.10)') 'Generalised Born:    ', A12DECOMPOSED_E % GB 
1665:                         WRITE(MUTUNIT,'(A,F20.10)') 'Surface energy:      ', A12DECOMPOSED_E % SURF 
1666:                         WRITE(MUTUNIT,'(A,F20.10)') 'Bond energy:         ', A12DECOMPOSED_E % BOND 
1667:                         WRITE(MUTUNIT,'(A,F20.10)') 'Angular term:        ', A12DECOMPOSED_E % ANGLE 
1668:                         WRITE(MUTUNIT,'(A,F20.10)') 'Dihedral term:       ', A12DECOMPOSED_E % DIHEDRAL 
1669:                         WRITE(MUTUNIT,'(A,F20.10)') 'vdW 1-4 term:        ', A12DECOMPOSED_E % VDW_14 
1670:                         WRITE(MUTUNIT,'(A,F20.10)') 'Electronic 1-4:      ', A12DECOMPOSED_E % ELEC_14 
1671:                         WRITE(MUTUNIT,'(A,F20.10)') 'Restraints:          ', A12DECOMPOSED_E % RESTRAINT 
1672:                         WRITE(MUTUNIT,'(A,F20.10)') 'Urey Bradley angle:  ', A12DECOMPOSED_E % ANGLE_UB 
1673:                         WRITE(MUTUNIT,'(A,F20.10)') 'Improper energy:     ', A12DECOMPOSED_E % IMP 
1674:                         WRITE(MUTUNIT,'(A,F20.10)') 'CMAP:                ', A12DECOMPOSED_E % CMAP 
1675:                      END IF1657:                      END IF
1676:                   END IF1658:                   END IF
1677:                ENDIF !<End of looong block of print statements1659:                ENDIF !<End of looong block of print statements
1678:                1660:                
1679: !     mp466>  writes structure and energetic data at regular increments1661: !     mp466>  writes structure and energetic data at regular increments
1680: !             to *plot and movie files for AMH potential1662: !             to *plot and movie files for AMH potential
1681: !     1663: !     
1682:                IF ((MOD(J1,NINT_AMH).EQ.0).AND.AMHT) THEN1664:                IF ((MOD(J1,NINT_AMH).EQ.0).AND.AMHT) THEN
1683:                   1665:                   
1684:                   GLY_COUNT = 01666:                   GLY_COUNT = 0
2396: 2378: 
2397:             !the mutational steps are accepted for better bonding properties, but the energies used 2379:             !the mutational steps are accepted for better bonding properties, but the energies used 
2398:             !otherwise refer to the best potential energies for one sequence2380:             !otherwise refer to the best potential energies for one sequence
2399:             IF (MUTATEDT.AND.((NMUTSTEP + MUTTESTSTEPS).EQ.J1)) THEN !kr366> here we check our mutational step            2381:             IF (MUTATEDT.AND.((NMUTSTEP + MUTTESTSTEPS).EQ.J1)) THEN !kr366> here we check our mutational step            
2400:                CALL MUTATION_E(EMUTN,BESTCOORDS,MUTENERGY,MUTTERMID) !get the energy for the new molecule2382:                CALL MUTATION_E(EMUTN,BESTCOORDS,MUTENERGY,MUTTERMID) !get the energy for the new molecule
2401:                WRITE(MYUNIT, '(2(A,F15.8))') ' mc> Energy before mutation: ',EMUTP,' | energy after mutation: ',EMUTN2383:                WRITE(MYUNIT, '(2(A,F15.8))') ' mc> Energy before mutation: ',EMUTP,' | energy after mutation: ',EMUTN
2402:                WRITE(MUTUNIT, '(2(A,F15.8))') 'Energy before mutation: ',EMUTP,' | energy after mutation: ',EMUTN2384:                WRITE(MUTUNIT, '(2(A,F15.8))') 'Energy before mutation: ',EMUTP,' | energy after mutation: ',EMUTN
2403:                IF (EMUTN.LT.EMUTP) THEN2385:                IF (EMUTN.LT.EMUTP) THEN
2404:                   RANDOM=0.0D02386:                   RANDOM=0.0D0
2405:                   ATEST=.TRUE.2387:                   ATEST=.TRUE.
2406:                ELSEIF ((MUTENERGY.EQ.2).OR.(MUTENERGY.EQ.3)) THEN !either one energy term or binding energy2388:                ELSE
2407:                   RANDOM=DPRAND()2389:                   RANDOM=DPRAND()
2408:                   DUMMY = EMUTN - EMUTP2390:                   DUMMY = EMUTN - EMUTP
2409:                   DUMMY=EXP(-DUMMY/MAX(MCTEMP,1.0D-100))2391:                   DUMMY=EXP(-DUMMY/MAX(MCTEMP,1.0D-100))
2410:                   IF (DUMMY.GT.RANDOM) THEN2392:                   IF (DUMMY.GT.RANDOM) THEN
2411:                      ATEST=.TRUE.2393:                      ATEST=.TRUE.
2412:                   ELSE2394:                   ELSE
2413:                      ATEST=.FALSE.2395:                      ATEST=.FALSE.
2414:                   ENDIF2396:                   ENDIF
2415:                ELSEIF (MUTENERGY.EQ.1) THEN  
2416:                   ATEST=.FALSE. !random choice 
2417:                ELSE !these are all distance based 
2418:                   RANDOM=DPRAND() 
2419:                   DUMMY = EMUTN/EMUTP 
2420:                   IF (DUMMY.GT.RANDOM) THEN 
2421:                      ATEST=.TRUE. 
2422:                   ELSE 
2423:                      ATEST=.FALSE. 
2424:                   ENDIF                   
2425:                ENDIF2397:                ENDIF
2426:                IF (ATEST) THEN2398:                IF (ATEST) THEN
2427:                   WRITE(MUTUNIT,'(A)') 'Mutation accepted'2399:                   WRITE(MUTUNIT,'(A)') 'Mutation accepted'
2428:                   WRITE(MYUNIT, '(A)') ' mc> Mutation accepted'2400:                   WRITE(MYUNIT, '(A)') ' mc> Mutation accepted'
2429:                   WRITE(MYUNIT,'(A)') ' mc> Sequence accepted: '2401:                   WRITE(MYUNIT,'(A)') ' mc> Sequence accepted: '
2430:                   CALL PRINT_CURRENT_SEQ()2402:                   CALL PRINT_CURRENT_SEQ()
2431:                ELSE2403:                ELSE
2432:                   CALL AMBER12_WRITE_RESTART(SCREENC,"before_reverse.mut",LEN("before_reverse.mut"))2404:                   CALL AMBER12_WRITE_RESTART(SCREENC,"before_reverse.mut",LEN("before_reverse.mut"))
2433:                   CALL SYSTEM("cp coords.prmtop before_mut.prmtop")2405:                   CALL SYSTEM("cp coords.prmtop before_mut.prmtop")
2434:                   CALL REVERSE_MUTATION(MUTATEDRES)2406:                   CALL REVERSE_MUTATION(MUTATEDRES)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0