hdiff output

r33469/amber_mutations.F90 2017-11-13 18:30:12.117626005 +0000 r33468/amber_mutations.F90 2017-11-13 18:30:12.341628969 +0000
611:      NSEQSTORED = NSEQSTORED + 1611:      NSEQSTORED = NSEQSTORED + 1
612:      MUTSEQ(NSEQSTORED, :) = AMBER12_RESNAME(:)612:      MUTSEQ(NSEQSTORED, :) = AMBER12_RESNAME(:)
613:      RETURN613:      RETURN
614:   END SUBROUTINE REVERSE_MUTATION614:   END SUBROUTINE REVERSE_MUTATION
615: 615: 
616:   !scoring function to judge how good mutation is616:   !scoring function to judge how good mutation is
617:   SUBROUTINE MUTATION_E(SCORE,COORDS,MODE,TERMID)617:   SUBROUTINE MUTATION_E(SCORE,COORDS,MODE,TERMID)
618:      DOUBLE PRECISION, INTENT(OUT) :: SCORE618:      DOUBLE PRECISION, INTENT(OUT) :: SCORE
619:      DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)619:      DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS)
620:      INTEGER, INTENT(IN) :: MODE, TERMID620:      INTEGER, INTENT(IN) :: MODE, TERMID
621:      DOUBLE PRECISION :: DPRAND, EREAL, NORMDIFF621:      DOUBLE PRECISION :: DPRAND, EREAL, GRADATOMS(3*NATOMS)
622:      DOUBLE PRECISION :: GRADATOMS(3*NATOMS), XSEP(3*NATOMS), XCA_1(3), XCA_2(3), DIFF12(3)622:      INTEGER :: ATOMID, PARMEDUNIT, GETUNIT, J1
623:      DOUBLE PRECISION, PARAMETER :: DSEP=1.5D2 
624:      INTEGER :: ATOMID, PARMEDUNIT, GETUNIT, J1, NDUMMY 
625:      TYPE(POT_ENE_REC_C) :: DECOMPOSED_E623:      TYPE(POT_ENE_REC_C) :: DECOMPOSED_E
626:      TYPE(AMBER12_ATOM) :: ATOMDATA(NATOMS) 
627:      CHARACTER(200) ENTRY_624:      CHARACTER(200) ENTRY_
628:      INTEGER , PARAMETER :: NWORDS=20625:      INTEGER , PARAMETER :: NWORDS=20
629:      CHARACTER(25) :: ENTRIES(NWORDS)=''626:      CHARACTER(25) :: ENTRIES(NWORDS)=''
630:      CHARACTER(LEN=6) :: J1_STRING627:      CHARACTER(LEN=6) :: J1_STRING
631:    628:    
632:      !random number as penalty function 
633:      IF (MODE.EQ.1) THEN629:      IF (MODE.EQ.1) THEN
634:         SCORE=DPRAND()630:         SCORE=DPRAND()
635:      !use contribution form AMBER energies 
636:      ELSE IF (MODE.EQ.2) THEN631:      ELSE IF (MODE.EQ.2) THEN
637:         CALL AMBER12_ENERGY_AND_GRADIENT(NATOMS, COORDS, EREAL, GRADATOMS, DECOMPOSED_E)632:         CALL AMBER12_ENERGY_AND_GRADIENT(NATOMS, COORDS, EREAL, GRADATOMS, DECOMPOSED_E)
638:         WRITE(MUTUNIT,'(A)') 'Energy decomposition'633:         WRITE(MUTUNIT,'(A)') 'Energy decomposition'
639:         WRITE(MUTUNIT,'(A,F20.10)') 'Total energy:        ', DECOMPOSED_E % TOTAL634:         WRITE(MUTUNIT,'(A,F20.10)') 'Total energy:        ', DECOMPOSED_E % TOTAL
640:         WRITE(MUTUNIT,'(A,F20.10)') 'Total van der Waals: ', DECOMPOSED_E % VDW_TOT635:         WRITE(MUTUNIT,'(A,F20.10)') 'Total van der Waals: ', DECOMPOSED_E % VDW_TOT
641:         WRITE(MUTUNIT,'(A,F20.10)') 'Total electronic:    ', DECOMPOSED_E % ELEC_TOT636:         WRITE(MUTUNIT,'(A,F20.10)') 'Total electronic:    ', DECOMPOSED_E % ELEC_TOT
642:         WRITE(MUTUNIT,'(A,F20.10)') 'Generalised Born:    ', DECOMPOSED_E % GB637:         WRITE(MUTUNIT,'(A,F20.10)') 'Generalised Born:    ', DECOMPOSED_E % GB
643:         WRITE(MUTUNIT,'(A,F20.10)') 'Surface energy:      ', DECOMPOSED_E % SURF638:         WRITE(MUTUNIT,'(A,F20.10)') 'Surface energy:      ', DECOMPOSED_E % SURF
644:         WRITE(MUTUNIT,'(A,F20.10)') 'Bond energy:         ', DECOMPOSED_E % BOND639:         WRITE(MUTUNIT,'(A,F20.10)') 'Bond energy:         ', DECOMPOSED_E % BOND
645:         WRITE(MUTUNIT,'(A,F20.10)') 'Angular term:        ', DECOMPOSED_E % ANGLE640:         WRITE(MUTUNIT,'(A,F20.10)') 'Angular term:        ', DECOMPOSED_E % ANGLE
672:            SCORE = DECOMPOSED_E % ELEC_14667:            SCORE = DECOMPOSED_E % ELEC_14
673:         ELSE IF (TERMID.EQ.10) THEN668:         ELSE IF (TERMID.EQ.10) THEN
674:            SCORE = DECOMPOSED_E % RESTRAINT669:            SCORE = DECOMPOSED_E % RESTRAINT
675:         ELSE IF (TERMID.EQ.11) THEN670:         ELSE IF (TERMID.EQ.11) THEN
676:            SCORE = DECOMPOSED_E % ANGLE_UB671:            SCORE = DECOMPOSED_E % ANGLE_UB
677:         ELSE IF (TERMID.EQ.12) THEN672:         ELSE IF (TERMID.EQ.12) THEN
678:            SCORE = DECOMPOSED_E % IMP673:            SCORE = DECOMPOSED_E % IMP
679:         ELSE IF (TERMID.EQ.13) THEN674:         ELSE IF (TERMID.EQ.13) THEN
680:            SCORE = DECOMPOSED_E % CMAP675:            SCORE = DECOMPOSED_E % CMAP
681:         ENDIF676:         ENDIF
682:      !approximate brinding interactions 
683:      ELSE IF (MODE.EQ.3) THEN677:      ELSE IF (MODE.EQ.3) THEN
684:         ATOMID = AMBER12_RESEND(TERMID) !get last atom in first group678:         ATOMID = AMBER12_RESEND(TERMID) !get last atom in first group
685:         !old way - very slow as it uses external script 
686:         !open new file to write parmed input679:         !open new file to write parmed input
687:         !PARMEDUNIT = GETUNIT()680:         PARMEDUNIT = GETUNIT()
688:         !OPEN(PARMEDUNIT,FILE='parmed_in',STATUS='NEW')681:         OPEN(PARMEDUNIT,FILE='parmed_in',STATUS='NEW')
689:         !WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') 'addExclusions @1'        682:         WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') 'addExclusions @1'        
690:         !DO J1=2,ATOMID683:         DO J1=2,ATOMID
691:         !   WRITE(J1_STRING,'(I6)') J1684:            WRITE(J1_STRING,'(I6)') J1
692:         !   WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING))685:            WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING))
693:         !ENDDO686:         ENDDO
694:         !WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ' @1'       687:         WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ' @1'       
695:         !DO J1=2,ATOMID-1688:         DO J1=2,ATOMID-1
696:         !   WRITE(J1_STRING,'(I6)') J1689:            WRITE(J1_STRING,'(I6)') J1
697:         !   WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING))690:            WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING))
698:         !ENDDO691:         ENDDO
699:         !WRITE(J1_STRING,'(I6)') ATOMID692:         WRITE(J1_STRING,'(I6)') ATOMID
700:         !WRITE(PARMEDUNIT,'(A)') ','//TRIM(ADJUSTL(J1_STRING))693:         WRITE(PARMEDUNIT,'(A)') ','//TRIM(ADJUSTL(J1_STRING))
701:         !WRITE(J1_STRING,'(I6)') ATOMID+1694:         WRITE(J1_STRING,'(I6)') ATOMID+1
702:         !WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') 'addExclusions @'//TRIM(ADJUSTL(J1_STRING))        695:         WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') 'addExclusions @'//TRIM(ADJUSTL(J1_STRING))        
703:         !DO J1=ATOMID+2,NATOMS696:         DO J1=ATOMID+2,NATOMS
704:         !   WRITE(J1_STRING,'(I6)') J1697:            WRITE(J1_STRING,'(I6)') J1
705:         !   WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING))698:            WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING))
706:         !ENDDO699:         ENDDO
707:         !WRITE(J1_STRING,'(I6)') ATOMID+1700:         WRITE(J1_STRING,'(I6)') ATOMID+1
708:         !WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ' @'//TRIM(ADJUSTL(J1_STRING))       701:         WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ' @'//TRIM(ADJUSTL(J1_STRING))       
709:         !DO J1=ATOMID+2,NATOMS-1702:         DO J1=ATOMID+2,NATOMS-1
710:         !   WRITE(J1_STRING,'(I6)') J1703:            WRITE(J1_STRING,'(I6)') J1
711:         !   WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING))704:            WRITE(PARMEDUNIT,'(A)',ADVANCE='NO') ','//TRIM(ADJUSTL(J1_STRING))
712:         !ENDDO705:         ENDDO
713:         !WRITE(J1_STRING,'(I6)') NATOMS-1706:         WRITE(J1_STRING,'(I6)') NATOMS-1
714:         !WRITE(PARMEDUNIT,'(A)') ','//TRIM(ADJUSTL(J1_STRING))707:         WRITE(PARMEDUNIT,'(A)') ','//TRIM(ADJUSTL(J1_STRING))
715:         !WRITE(PARMEDUNIT,'(A)') 'loadRestrt current.inpcrd'708:         WRITE(PARMEDUNIT,'(A)') 'loadRestrt current.inpcrd'
716:         !WRITE(J1_STRING,'(I6)') AMBERMUTIGB709:         WRITE(J1_STRING,'(I6)') AMBERMUTIGB
717:         !WRITE(PARMEDUNIT,'(A)') 'energy cutoff 15.0 igb '//TRIM(ADJUSTL(J1_STRING))//' saltcon 0.1'710:         WRITE(PARMEDUNIT,'(A)') 'energy cutoff 15.0 igb '//TRIM(ADJUSTL(J1_STRING))//' saltcon 0.1'
718:         !WRITE(PARMEDUNIT,'(A)') 'quit'711:         WRITE(PARMEDUNIT,'(A)') 'quit'
719:         !CLOSE(PARMEDUNIT)712:         CLOSE(PARMEDUNIT)
720:         !CALL AMBER12_WRITE_RESTART(COORDS, 'current.inpcrd',LEN('current.inpcrd'))713:         CALL AMBER12_WRITE_RESTART(COORDS, 'current.inpcrd',LEN('current.inpcrd'))
721:         !!create new topology without interactions and calculate energy714:         !create new topology without interactions and calculate energy
722:         !CALL SYSTEM('parmed.py -n coords.prmtop parmed_in > parmed_out')715:         CALL SYSTEM('parmed.py -n coords.prmtop parmed_in > parmed_out')
723:         !OPEN(PARMEDUNIT,FILE='parmed_out',STATUS='OLD')716:         OPEN(PARMEDUNIT,FILE='parmed_out',STATUS='OLD')
724:         !DO717:         DO
725:         !  ENTRIES(:)=''718:           ENTRIES(:)=''
726:         !  READ(PARMEDUNIT,'(A)',END=588) ENTRY_719:           READ(PARMEDUNIT,'(A)',END=588) ENTRY_
727:         !  CALL READ_LINE(ENTRY_,NWORDS,ENTRIES)720:           CALL READ_LINE(ENTRY_,NWORDS,ENTRIES)
728:         !  IF (ENTRIES(1).EQ.'TOTAL') THEN721:           IF (ENTRIES(1).EQ.'TOTAL') THEN
729:         !     READ(ENTRIES(3),'(F20.10)') SCORE722:              READ(ENTRIES(3),'(F20.10)') SCORE
730:         !     GOTO 588723:              GOTO 588
731:         !  ENDIF724:           ENDIF
732:         !ENDDO725:         ENDDO
733: !588     CONTINUE726: 588     CONTINUE
734:         !CLOSE(PARMEDUNIT)727:         CLOSE(PARMEDUNIT)
735:         !CALL SYSTEM('rm current.inpcrd parmed_in parmed_out')728:         CALL SYSTEM('rm current.inpcrd parmed_in parmed_out')      
736:  
737:  
738:         !new way - simply move first and second half apart and compute amber energy 
739:         CALL AMBER12_GET_ATOMDATA(ATOMDATA,NATOMS) 
740:         !get the positions of the C alpha atoms and then assign the centre of these coordinates 
741:         NDUMMY = 0 
742:         XCA_1(:) = 0.0D0 
743:         XCA_2(:) = 0.0D0  
744:         DO J1=1,NATOMS 
745:             IF (ATOMDATA(J1)%NAME.EQ.'C') THEN 
746:                 NDUMMY = NDUMMY + 1 
747:                 IF (J1.LE.ATOMID) THEN 
748:                     XCA_1(1) = XCA_1(1) + COORDS(3*J1-2) 
749:                     XCA_1(2) = XCA_1(2) + COORDS(3*J1-1) 
750:                     XCA_1(3) = XCA_1(3) + COORDS(3*J1) 
751:                     IF (NDUMMY.EQ.TERMID) THEN 
752:                         XCA_1(1) = XCA_1(1)/TERMID 
753:                         XCA_1(1) = XCA_1(1)/TERMID 
754:                         XCA_1(1) = XCA_1(1)/TERMID 
755:                     ENDIF 
756:                 ELSE 
757:                     XCA_2(1) = XCA_2(1) + COORDS(3*J1-2) 
758:                     XCA_2(2) = XCA_2(2) + COORDS(3*J1-1) 
759:                     XCA_2(3) = XCA_2(3) + COORDS(3*J1) 
760:                 ENDIF 
761:             ENDIF 
762:             IF (NDUMMY.EQ.NRESIDUES) GOTO 760 
763:         END DO 
764: 760     CONTINUE 
765:         XCA_1(2) = XCA_1(2)/(NRESIDUES-TERMID) 
766:         XCA_1(2) = XCA_1(2)/(NRESIDUES-TERMID) 
767:         XCA_1(2) = XCA_1(2)/(NRESIDUES-TERMID) 
768:         !now define the vector between the two centres 
769:         DIFF12(1) = XCA_2(1) - XCA_1(1) 
770:         DIFF12(2) = XCA_2(2) - XCA_1(2) 
771:         DIFF12(3) = XCA_2(3) - XCA_1(3) 
772:         NORMDIFF = 1.0D0/SQRT(DIFF12(1)**2 + DIFF12(2)**2 + DIFF12(3)**2) 
773:         DIFF12(1) = NORMDIFF * DIFF12(1) 
774:         DIFF12(2) = NORMDIFF * DIFF12(2) 
775:         DIFF12(3) = NORMDIFF * DIFF12(3) 
776:         !now we have the normalised vector pointing along the Calpha centres 
777:         XSEP(:) = COORDS(:) 
778:         DO J1=(ATOMID+1),NATOMS 
779:            XSEP(3*J1-2) = COORDS(3*J1-2) + DSEP * DIFF12(1) 
780:            XSEP(3*J1-1) = COORDS(3*J1-1) + DSEP * DIFF12(1) 
781:            XSEP(3*J1) = COORDS(3*J1) + DSEP * DIFF12(1) 
782:        END DO 
783:        CALL AMBER12_ENERGY_AND_GRADIENT(NATOMS, XSEP, SCORE, GRADATOMS, DECOMPOSED_E) 
784:  
785:      ELSE IF (MODE.EQ.4) THEN729:      ELSE IF (MODE.EQ.4) THEN
786:          IF (.NOT.(ALLOCATED(CA_REFERENCE))) THEN730:          IF (.NOT.(ALLOCATED(CA_REFERENCE))) THEN
787:             ALLOCATE(CA_REFERENCE(3*NRESIDUES))731:             ALLOCATE(CA_REFERENCE(3*NRESIDUES))
788:             CALL CREATE_CA_REF(TERMID)732:             CALL CREATE_CA_REF(TERMID)
789:          END IF733:          END IF
790:          CALL CALPHA_RMSD(SCORE,COORDS)734:          CALL CALPHA_RMSD(SCORE,COORDS)
791:      !helical optimisation for alpha helix735:      !helical optimisation for alpha helix
792:      ELSE IF (MODE.EQ.5) THEN736:      ELSE IF (MODE.EQ.5) THEN
793:          CALL HELICAL_DSSP(SCORE,COORDS,1)737:          CALL HELICAL_DSSP(SCORE,COORDS,1)
794:      !helical optimisation for 3-10 helix738:      !helical optimisation for 3-10 helix


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0