diff --git a/contrib/pdb_tools/GNUmakefile b/contrib/pdb_tools/GNUmakefile new file mode 100644 index 0000000000..5b56052ade --- /dev/null +++ b/contrib/pdb_tools/GNUmakefile @@ -0,0 +1,7 @@ +OBJ_OPTIMIZE=pdb_scale.o + +all: $(OBJ_OPTIMIZE) + $(FC) -o pdb_scale $(OBJ_OPTIMIZE) + +.F.o: + $(FC) -c $< diff --git a/contrib/pdb_tools/README.md b/contrib/pdb_tools/README.md new file mode 100644 index 0000000000..8cee4121fa --- /dev/null +++ b/contrib/pdb_tools/README.md @@ -0,0 +1,33 @@ +Tools to manipulate PDB files +============================= + +General PDB files need tools to apply transformations to them. As this need is +quite common there a number of tools out there already and where applicable +you should use them. Examples of such tools are: + +* convpdb.pl +* pdb4amber +* charmmlipid2amber.py +* pdbfixer + +Unfortunately some of these tools try to be too helpful. For example `convpdb.pl` +will not only scale the coordinates in a PDB file if you ask for it, but it will +also drop all TER records, rename HOH to TIP3, drop the CRYST1 record, and drop +the chemical symbol at the end of the line. This might leave you with a useless +result. + +The tools in this directory do exactly what it says on the tin, nothing less +and nothing more. The tools provided here are: + +* pdb_amber2nwchem + * Convert an Amber PDB file into an NWChem PDB file by renaming atoms and + residues. +* pdb_large + * Interconvert a PDB file between the regular PDB file format and the + NWChem PDB file format for large molecular systems. In essence the + large PDB format allows 6 digits for the residue number instead of 4. +* pdb_scale.F + * Scale the atomic coordinates and the box size with a specified factor. +* pdb_supercell.py + * Build a supercell of a given PDB file by doubling the initial molecular + system in the requested direction. diff --git a/contrib/pdb_tools/pdb_amber2nwchem/GNUmakefile b/contrib/pdb_tools/pdb_amber2nwchem/GNUmakefile new file mode 100644 index 0000000000..5229bea7fa --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/GNUmakefile @@ -0,0 +1,14 @@ +OBJ_OPTIMIZE=inp.o \ + pdb_amber2nwchem.o \ + pdb_readcommandlinearguments.o \ + pre_atnum.o \ + pre_namiup.o \ + pre_protein.o \ + util_print.o + + +all: $(OBJ_OPTIMIZE) + $(FC) -o pdb_amber2nwchem $(OBJ_OPTIMIZE) + +.F.o: + $(FC) -c $< diff --git a/contrib/pdb_tools/pdb_amber2nwchem/README.md b/contrib/pdb_tools/pdb_amber2nwchem/README.md new file mode 100644 index 0000000000..b226cd7a72 --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/README.md @@ -0,0 +1,38 @@ +# README PDB_AMBER2NWCHEM + +`pdb_amber2nwchem` is a tool to convert PDB files from the Amber +dialect into the NWChem dialect. PDB files have standardized atom +names for many residues of biological relevance. However, most PDB +files do not contain Hydrogen atoms and different conventions have +been adopted by different code communities. This tool converts +a PDB file in the Amber convention to one following the NWChem +conventions, so that NWChem can read it successfully. + +Apart from atom names the tool also converts some residue names: +- "WAT" --> "HOH" +- "Na+" --> "Na" +- "K+" --> "K" +- "Cl-" --> "Cl" + +to name a few examples. + +# References + +- [1] "IUPAC-IUB Commission on Biochemical Nomenclature. + Abbreviations and Symbols for the Description of the + Conformation of Polypeptide Chains. Tentative Rules (1969)". + Biochemistry (1970) 9, 3471-3479. + DOI: [10.1021/bi00820a001](https://doi.org/10.1021/bi00820a001) + +- [2] Charles Hoogstraten, "Correlation of hydrogen atom naming + systems, including diastereotopic protons." Web-page: + http://www.bmrb.wisc.edu/ref_info/atom_nom.tbl + [accessed Jan 6, 2018] + +- [3] IUPAC-IUBMB Joint Commission on Biochemical Nomenclature + and Nomenclature Commission of IUBMB, "Biochemical + Nomenclature and Related Documents (White Book)." 2nd Edition, + Portland Press, 1992, pp. 39-69. ISBN: 1-85578-005-4. + Part: "Amino acids, peptides and proteins", Section 3AA-1 + "Names of Common alpha-Amino Acids". + DOI: [10.1351/pac198456050595](https://doi.org/10.1351/pac198456050595) diff --git a/contrib/pdb_tools/pdb_amber2nwchem/inp.F b/contrib/pdb_tools/pdb_amber2nwchem/inp.F new file mode 100644 index 0000000000..b59d52620c --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/inp.F @@ -0,0 +1,61 @@ + block data inp_data +C$Id$ + implicit none +#include "inpP.fh" +c + data iread /5/ + data iwrite /6/ + data jrec /-1/ + data jump /0/ + data oswit /.false. / + data nerr /999/ + data nline /0/ + data noline /0/ + data ierrpos /-1/ + data errmsg /' '/ + data input_line /0/ + data xblnk /' '/ +#if defined(SGITFP) || defined(SGI_N32) || defined(CRAY) || defined(HPUX) || defined(WIN32) || defined(PSCALE) || defined(__FLANG) || ( __GNUC__ >= 4) + data xtab /' '/ ! Tab ... no backslash necessary +#elif (defined(LINUX) || defined(MACX)) && !defined(PGLINUX) && !defined(XLFLINUX) &&!( __GNUC__ >= 4) + data xtab /9/ ! Tab ... g77 has trouble with escape sequence +#else + data xtab /'\ '/ ! Tab ... note backslash for cpp +#endif + data xsplit/';'/ + data xcomm /'#'/ + data xback /'\\'/ ! Backslash ... note backslash for cpp + data xquote/'"'/ +c + data save_level /0/ ! Depth of nesting in save/restore + data include_level /0/ +c + end +c +C> \brief Returns the length of the contents of a string +C> +C> In Fortran the length of a string is determined by the amount of +C> memory allocated for it. Often one just wants to know how long the +C> value of the string is without any trailing spaces. Which is what +C> this function returns. +C> +C> \return The length of the value of a string without any trailing +C> spaces +c + integer function inp_strlen(a) + implicit none +#include "inpP.fh" + character*(*) a !< [Input] The string + integer i + logical ois_ws + intrinsic len + character*1 xtest + ois_ws(xtest) = (xtest.eq.xblnk .or. xtest.eq.xtab) +c + do i = len(a),1,-1 + if (.not. ois_ws(a(i:i))) goto 10 + enddo +c + 10 inp_strlen = i +c + end diff --git a/contrib/pdb_tools/pdb_amber2nwchem/inpP.fh b/contrib/pdb_tools/pdb_amber2nwchem/inpP.fh new file mode 100644 index 0000000000..9948bd862a --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/inpP.fh @@ -0,0 +1,76 @@ +c +c Private header file for free format input routines +c +** implicit real*8 (a-h,p-w), integer (i-n), logical (o) +** implicit character *8 (z),character *1 (x) +** implicit character *4 (y) +c + integer max_width ! Maximum input width of an input line +C$Id$ + integer max_field ! Maximum no. of fields in an input line + parameter (max_width = 1024) + parameter (max_field = max_width/2 + 1) + character*1024 ja, ia ! Input buffers ... MUST match max_width + character*1024 tmp ! Same size work space + character*80 errmsg ! Error message + character*1 xcomm ! Comment character + character*1 xsplit ! Character to split physical input lines + character*1 xback ! Backslash for concatenation and quoting + character*1 xquote ! Quotation marks for strings + character*1 xblnk ! Space + character*1 xtab ! Tab + common /inp_workc / ia, ja, tmp, errmsg, + $ xcomm, xsplit, xback, xquote, xblnk, xtab +c + integer iread, iwrite + integer jrec ! No. of current field + integer jump ! No. of fields in current line + integer istrt ! Start of fields + integer inumb ! Length of fields + integer nend ! End of fields + integer iwidth ! Length of current logical input line + integer nline ! Current logical line inside physical line + integer noline ! No. of logical lines inside physical line + integer input_line ! No. of current physical input line + integer nerr ! ? + logical oswit ! True if EOF has beeen detected + integer ierrpos ! Input char position where error was detected + integer nstart +c + common /inp_work/ jrec, jump, + $ istrt(max_field), inumb(max_field),iwidth,nend(max_field), + $ nstart(max_field),nline,noline,nerr,oswit,ierrpos, input_line + common/inp_file/iread,iwrite +c + integer max_include_level + parameter (max_include_level = 3) + integer include_level + common/inp_include_info/include_level + character*255 include_file_name + common/inp_include_cinfo/include_file_name(max_include_level) +c +c This drivel below is to save state so that by doing the following +c +c call inp_save_state +c +c call inp_init +c +c call inp_restore_state +c +c It is possible to go off and read another file without +c destroying the state of the current file. Note that this +c only works nested three deep (i.e., 4 files concurrently open) +c +c Note that all of these commons must match in size the +c ones above and that inp_save/restore_state must also be +c consistent +c + character*1 csave + integer isave + integer fsave + integer save_level + common/inp_csave/csave(3*max_width+80+6,3) + common/inp_isave/isave(9+4*max_field,3) + common/inp_fsave/fsave(2,3) + common/inp_save_data/save_level +c diff --git a/contrib/pdb_tools/pdb_amber2nwchem/pdb_amber2nwchem.F b/contrib/pdb_tools/pdb_amber2nwchem/pdb_amber2nwchem.F new file mode 100644 index 0000000000..572b43dfb6 --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/pdb_amber2nwchem.F @@ -0,0 +1,101 @@ + program pdb_amber2nwchem + implicit none +#include "pre_common.fh" + character(len=1024) :: inputfile + character(len=1024) :: outputfile + character(len=100) :: line + character(len=12) :: l1 ! characters 1-12 + character(len=6) :: catm(3,1) ! characters 13-16 + character(len=1) :: l2 ! characters 17 + character(len=5) :: cseq(2,1) ! characters 18-20 + character(len=100) :: l3 ! characters 21-72 + integer, external :: inp_strlen + logical, external :: pre_namiup + integer :: lfninp + integer :: lfnout + integer :: iline + integer :: lseq(6,1) + integer :: latm(5,1) + logical :: ll + ffield = "amber" + call pdb_readcommandlinearguments(inputfile,outputfile,ffield) + lfninp = 5 + lfnout = 6 + lseq = 1 + latm = 1 + if (inputfile.ne."-") then + close(lfninp) + open(unit=lfninp,file=inputfile,err=100,form="formatted", + + status="old") + endif + 100 continue + if (outputfile.ne."-") then + close(lfnout) + open(unit=lfnout,file=outputfile,err=200,form="formatted") + endif + 200 continue + iline = 0 + line = "" + do while (line(1:3).ne."END") + catm = " " + cseq = " " + iline = iline + 1 + read(lfninp,'(a)')line + select case (line(1:4)) + case ("ATOM ", "HETA") + read(line,1002,end=300,err=300)l1,catm(1,1)(1:4),l2, + + cseq(1,1)(1:3),l3 + ll = pre_namiup(0,lseq,cseq,1,1,latm,catm,1,1) + call namseq(catm,cseq) + write(lfnout,1002)l1,catm(2,1)(1:4),l2,cseq(2,1)(1:3), + + l3(1:inp_strlen(l3)) + case default + write(lfnout,'(a)')line(1:inp_strlen(line)) + end select + enddo + 300 continue + close(lfninp) + close(lfnout) + 1002 format(a12,a4,a1,a3,a) + end +! +!----------------------------------------------------------------------- +! + subroutine namseq(catm,cseq) + implicit none + character(len=6), intent(inout) :: catm(3,1) + character(len=5), intent(inout) :: cseq(2,1) +! + cseq(2,1)(1:3) = cseq(1,1)(1:3) + if (cseq(2,1)(1:3).eq."WAT") then + cseq(2,1)(1:3) = "HOH" + if (catm(2,1)(1:4).eq." H1 ") then + catm(2,1)(1:4) = "2H " + endif + if (catm(2,1)(1:4).eq." H2 ") then + catm(2,1)(1:4) = "3H " + endif + if (catm(2,1)(1:4).eq." H3 ") then + catm(2,1)(1:4) = "2H " + endif + endif + if (cseq(2,1)(1:3).eq."Cl-") then + cseq(2,1)(1:3) = " Cl" + if (catm(2,1)(1:4).eq."Cl- ") then + catm(2,1)(1:4) = "Cl " + endif + endif + if (cseq(2,1)(1:3).eq."Na+") then + cseq(2,1)(1:3) = " Na" + if (catm(2,1)(1:4).eq."Na+ ") then + catm(2,1)(1:4) = "Na " + endif + endif + if (cseq(2,1)(1:3).eq." K+") then + cseq(2,1)(1:3) = " K " + if (catm(2,1)(1:4).eq." K+ ") then + catm(2,1)(1:4) = " K " + endif + endif + return + end diff --git a/contrib/pdb_tools/pdb_amber2nwchem/pdb_readcommandlinearguments.F b/contrib/pdb_tools/pdb_amber2nwchem/pdb_readcommandlinearguments.F new file mode 100644 index 0000000000..6ed973d5aa --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/pdb_readcommandlinearguments.F @@ -0,0 +1,64 @@ + subroutine pdb_readcommandlinearguments(pdb_input,pdb_output, + + pdb_dialect) + implicit none + character(len=*), intent(out) :: pdb_input + character(len=*), intent(out) :: pdb_output + character(len=*), intent(out) :: pdb_dialect +! + integer :: narg ! the number of command line arguments + integer :: iarg ! the current command line argument +! + character(len=8) :: arg +! + pdb_input = "-" + pdb_output = "-" + narg = command_argument_count() + select case (narg) + case (2, 4, 6) + iarg=1 + do while (iarg.lt.narg) + call get_command_argument(iarg,arg) + select case (arg) + case ("-i") + iarg=iarg+1 + call get_command_argument(iarg,pdb_input) + case ("-o") + iarg=iarg+1 + call get_command_argument(iarg,pdb_output) + case ("-f") + iarg=iarg+1 + call get_command_argument(iarg,arg) + if (arg(1:5).eq."amber") then + pdb_dialect = "amber" + else if (arg(1:6).eq."charmm") then + pdb_dialect = "charmm" + else + write(*,*)"Error: unsupported forcefield: ",arg + write(*,*)'Supported forcefields: "amber" or + +"charmm"' + write(*,*) + go to 10 + endif + case default + go to 10 + end select + iarg=iarg+1 + enddo + case default + go to 10 + end select + return +! + 10 write(*,*)'Usage: pdb_amber2nwchem -i inputfile -o outputfile' + write(*,*)' [-f forcefield]' + write(*,*) + write(*,*)'inputfile : the name of the Amber PDB input file.' + write(*,*) + write(*,*)'outputfile: the name of the NWChem PDB output + +file.' + write(*,*) + write(*,*)'forcefield: the name of the forcefield (default + +"amber").' + stop 10 +! + end diff --git a/contrib/pdb_tools/pdb_amber2nwchem/pre_atnum.F b/contrib/pdb_tools/pdb_amber2nwchem/pre_atnum.F new file mode 100644 index 0000000000..5580fd10ca --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/pre_atnum.F @@ -0,0 +1,51 @@ + integer function pre_atnum(element) +c +c $Id$ +c + implicit none + character*2 element +c + character*2 aname(0:105),bname(0:105) +c + data aname / ' ', + + ' H','He','Li','Be',' B',' C',' N',' O',' F','Ne', + + 'Na','Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca', + + 'Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn', + + 'Ga','Ge','As','Se','Br',' R','Rb','Sr',' Y','Zr', + + 'Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn', + + 'Sb','Te',' I','Xe','Cs','Ba','La','Ce','Pr','Nd', + + 'Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb', + + 'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg', + + 'Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th', + + 'Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm', + + 'Md','No','Lr','Rf','Ha' / + data bname / ' ', + + ' H','HE','LI','BE',' B',' C',' N',' O',' F','NE', + + 'NA','MG','AL','SI',' P',' S','CL','AR',' K','CA', + + 'SC','TI',' V','CR','MN','FE','CO','NI','CU','ZN', + + 'GA','GE','AS','SE','BR',' R','RB','SR',' Y','ZR', + + 'NB','MO','TC','RU','RH','PD','AG','CD','IN','SN', + + 'SB','TE',' I','XE','CS','BA','LA','CE','PR','ND', + + 'PM','SM','EU','GD','TB','DY','HO','ER','TM','YB', + + 'LU','HF','TA',' W','RE','OS','IR','PT','AU','HG', + + 'TL','PB','BI','PO','AT','RN','FR','RA','AC','TH', + + 'PA',' U','Np','PU','AM','CM','BK','CF','ES','FM', + + 'MD','No','LR','RF','HA' / +c + integer i +c + do 1 i=1,105 + if(element.eq.aname(i).or.element.eq.bname(i)) then + pre_atnum=i + return + endif + 1 continue +c + pre_atnum=0 + if(element.eq.'1H') pre_atnum=1 + if(element.eq.'2H') pre_atnum=1 + if(element.eq.'3H') pre_atnum=1 + if(element.eq.'4H') pre_atnum=1 +c + return + end diff --git a/contrib/pdb_tools/pdb_amber2nwchem/pre_common.fh b/contrib/pdb_tools/pdb_amber2nwchem/pre_common.fh new file mode 100644 index 0000000000..a0e1389259 --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/pre_common.fh @@ -0,0 +1,50 @@ +c $Id$ +c +c general parameters +c + real*8 zero,tiny,small + real*8 fourth,third,half,one,two,three,four,six,twelve + parameter(zero=0.0d0) + parameter(tiny=1.0d-12) + parameter(small=1.0d-6) + parameter(one=1.0d0) + parameter(fourth=2.5d-1) + parameter(third=3.333333333333d-1) + parameter(half=5.0d-1) + parameter(two=2.0d0) + parameter(three=3.0d0) + parameter(four=4.0d0) + parameter(six=6.0d0) + parameter(twelve=1.2d1) +c +c dimensions +c +c mdirpar : maximum number of parameter directories +c mfilpar : maximum number of parameter files per parameter directory +c + integer mdirpar,mfilpar + parameter(mdirpar=6) + parameter(mfilpar=4) +c +c ndirpar : number of parameter directories +c nfilpar(i) : number of parameter files in i-th directory +c +c ffield : name of the force field +c dirpar(i) : name of i-th parameter directory +c filpar(i,j) : name of j-th parameter file in i-th directory +c + integer ndirpar,nfilpar(mdirpar) +c + common/pre_dim/ndirpar,nfilpar +c + character*80 ffield + character*255 dirpar(mdirpar),filpar(mdirpar,mfilpar) + character*255 filmst +c + common/pre_chr/ffield,dirpar,filpar,filmst +c + integer lfnmst + common/pre_lfn/lfnmst +c + logical lqhopmem + common/pre_log/lqhopmem diff --git a/contrib/pdb_tools/pdb_amber2nwchem/pre_namiup.F b/contrib/pdb_tools/pdb_amber2nwchem/pre_namiup.F new file mode 100644 index 0000000000..ddb6114eba --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/pre_namiup.F @@ -0,0 +1,2077 @@ + logical function pre_namiup(lfnout, + + lseq,cseq,mseq,nseq,latm,catm,matm,natm) +c +c $Id$ +c +c Function to comply with IUPAC atom naming [1-3]. +c +c in : lfnout = output file logical file number +c kse +c mseq = dimension of the sequence list +c nseq = length of the sequence list +c +c out : lseq(1,*) = segment numbers +c lseq(2,*) = number of atoms in segment +c lseq(3,*) = index to list of unique segments +c cseq = segment names +c nseq = actual length of the sequence list +c +c References +c +c [1] "IUPAC-IUB Commission on Biochemical Nomenclature. +c Abbreviations and Symbols for the Description of the +c Conformation of Polypeptide Chains. Tentative Rules (1969)". +c Biochemistry (1970) 9, 3471-3479. DOI: 10.1021/bi00820a001 +c +c [2] Charles Hoogstraten, "Correlation of hydrogen atom naming +c systems, including diastereotopic protons." Web-page: +c http://www.bmrb.wisc.edu/ref_info/atom_nom.tbl +c [accessed Jan 6, 2018] +c +c [3] IUPAC-IUBMB Joint Commission on Biochemical Nomenclature +c and Nomenclature Commission of IUBMB, "Biochemical +c Nomenclature and Related Documents (White Book)." 2nd Edition, +c Portland Press, 1992, pp. 39-69. ISBN: 1-85578-005-4. +c Part: "Amino acids, peptides and proteins", Section 3AA-1 +c "Names of Common alpha-Amino Acids". +c DOI: 10.1351/pac198456050595 +c + implicit none +c +#include "util.fh" +#include "pre_common.fh" +c + integer pre_atnum + external pre_atnum + logical pre_protein + external pre_protein +c + integer lfnout + integer mseq,matm + integer nseq,natm + integer lseq(6,mseq),latm(5,matm) + character*6 catm(3,matm) + character*5 cseq(2,mseq) + character*2 elemnt + character*3 cres ! residue name +c + integer i + logical flag1 +c + if(util_print('sequence',print_debug)) then + write(lfnout,2000) + 2000 format(/,'pre_namiup ') + endif +c + flag1=.false. + do 1 i=1,natm + if(pre_protein(cseq(1,latm(1,i)))) then + if(catm(1,i)(1:2).eq.'1H') flag1=.true. + endif + catm(2,i)=catm(1,i) + 1 continue +c +c IUPAC hydrogen naming +c + if(flag1) then + do 2 i=1,natm + if(pre_protein(cseq(1,latm(1,i)))) then + if(catm(2,i)(1:2).eq.'3H') catm(2,i)(1:2)='4H' + if(catm(2,i)(1:2).eq.'2H') catm(2,i)(1:2)='3H' + if(catm(2,i)(1:2).eq.'1H') catm(2,i)(1:2)='2H' + endif + 2 continue + if(util_print('sequence',print_high)) then + write(lfnout,2001) + 2001 format(/,'Hydrogen atom prefix corrected ') + endif + endif +c + 2002 format('Renamed ',a3,i3.3,x,a) +c +c ALA HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 4HB +c + cres='ALA' + do 3 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'ALA') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='4HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 4HB' + endif + endif + endif + 3 continue +c +c ARG HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c HG1 -> 2HG +c HG2 -> 3HG +c HG3 -> 2HG +c HD1 -> 2HD +c HD2 -> 3HD +c HD3 -> 2HD +c HH11 -> 2HH1 +c HH12 -> 3HH1 +c HH21 -> 2HH2 +c HH22 -> 3HH2 +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. Similary for HG1, HG2 and HD1, HD2. +c + cres='ARG' + do 4 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'ARG') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HG1') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG1 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG2') then + catm(2,i)(1:4)='3HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG2 to 3HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG3') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG3 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.' HD1') then + catm(2,i)(1:4)='2HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD1 to 2HD' + endif + endif + if(catm(2,i)(1:4).eq.' HD2') then + catm(2,i)(1:4)='3HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD2 to 3HD' + endif + endif + if(catm(2,i)(1:4).eq.' HD3') then + catm(2,i)(1:4)='2HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD3 to 2HD' + endif + endif + if(catm(2,i)(1:4).eq.'HH11') then + catm(2,i)(1:4)='2HH1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HH11 to 2HH1' + endif + endif + if(catm(2,i)(1:4).eq.'HH12') then + catm(2,i)(1:4)='3HH1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HH12 to 3HH1' + endif + endif + if(catm(2,i)(1:4).eq.'HH21') then + catm(2,i)(1:4)='2HH2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HH21 to 2HH2' + endif + endif + if(catm(2,i)(1:4).eq.'HH22') then + catm(2,i)(1:4)='3HH2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HH22 to 3HH2' + endif + endif + endif + 4 continue +c +c ASN HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c HD21 -> 2HD2 +c HD22 -> 3HD2 +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. +c + cres='ASN' + do 5 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'ASN') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.'HD21') then + catm(2,i)(1:4)='2HD2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD21 to 2HD2' + endif + endif + if(catm(2,i)(1:4).eq.'HD22') then + catm(2,i)(1:4)='3HD2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD22 to 3HD2' + endif + endif + endif + 5 continue +c +c ASP HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. +c + cres='ASP' + do 6 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'ASP') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + endif + 6 continue +c +c CYS HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. +c + cres='CYS' + do 7 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'CYS'.or. + + cseq(1,latm(1,i))(1:3).eq.'CYX') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + endif + 7 continue +c +c GLN HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c HG1 -> 2HG +c HG2 -> 3HG +c HG3 -> 2HG +c HE21 -> 2HE2 +c HE22 -> 3HE2 +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. Likewise for the HG atoms. +c + cres='GLN' + do 8 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'GLN') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HG1') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG1 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG2') then + catm(2,i)(1:4)='3HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG2 to 3HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG3') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG3 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.'HE21') then + catm(2,i)(1:4)='2HE2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HE21 to 2HE2' + endif + endif + if(catm(2,i)(1:4).eq.'HE22') then + catm(2,i)(1:4)='3HE2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HE22 to 3HE2' + endif + endif + endif + 8 continue +c +c GLU HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c HG1 -> 2HG +c HG2 -> 3HG +c HG3 -> 2HG +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. Likewise for the HG atoms. +c + cres='GLU' + do 9 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'GLU') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HG1') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG1 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG2') then + catm(2,i)(1:4)='3HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG2 to 3HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG3') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG3 to 2HG' + endif + endif + endif + 9 continue +c +c GLY HN -> H +c HA1 -> 2HA +c HA2 -> 3HA +c HA3 -> 2HA +c +c Some PDB files list HA1, HA2, some list HA2, HA3, ultimately both +c cases need to map to 2HA, 3HA. Hence the funky mapping of HA1, +c HA2, and HA3. +c + cres='GLY' + do 10 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'GLY') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HA1') then + catm(2,i)(1:4)='2HA ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HA1 to 2HA' + endif + endif + if(catm(2,i)(1:4).eq.' HA2') then + catm(2,i)(1:4)='3HA ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HA2 to 3HA' + endif + endif + if(catm(2,i)(1:4).eq.' HA3') then + catm(2,i)(1:4)='2HA ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HA3 to 2HA' + endif + endif + endif + 10 continue +c +c HIS HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. +c + cres='HIS' + do 11 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'HIS'.or. + + cseq(1,latm(1,i))(1:3).eq.'HID'.or. ! amber + + cseq(1,latm(1,i))(1:3).eq.'HIE'.or. + + cseq(1,latm(1,i))(1:3).eq.'HIP'.or. + + cseq(1,latm(1,i))(1:3).eq.'HSD'.or. ! charmm + + cseq(1,latm(1,i))(1:3).eq.'HSE'.or. + + cseq(1,latm(1,i))(1:3).eq.'HSP') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + endif + 11 continue +c +c ILE HN -> H +c HG21 -> 2HG2 +c HG22 -> 3HG2 +c HG23 -> 4HG2 +c HG11 -> 2HG1 +c HG12 -> 3HG1 +c HG13 -> 2HG1 +c CD1 -> CD +c HD1 -> 2HD +c HD2 -> 3HD +c HD3 -> 4HD +c HD11 -> 2HD +c HD12 -> 3HD +c HD13 -> 4HD +c +c Some PDB files list HG11, HG12, some list HG12, HG13, ultimately +c both cases need to map to 2HG1, 3HG1. Hence the funky mapping of +c HG11, HG12, and HG13. +c Also some PDB files list HD1, HD2, and HD3, some instead list +c HD11, HD12, and HD13. +c + cres='ILE' + do 12 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'ILE') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.'HG21') then + catm(2,i)(1:4)='2HG2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG21 to 2HG2' + endif + endif + if(catm(2,i)(1:4).eq.'HG22') then + catm(2,i)(1:4)='3HG2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG22 to 3HG2' + endif + endif + if(catm(2,i)(1:4).eq.'HG23') then + catm(2,i)(1:4)='4HG2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG22 to 4HG2' + endif + endif + if(catm(2,i)(1:4).eq.'HG11') then + catm(2,i)(1:4)='2HG1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG11 to 2HG1' + endif + endif + if(catm(2,i)(1:4).eq.'HG12') then + catm(2,i)(1:4)='3HG1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG12 to 3HG1' + endif + endif + if(catm(2,i)(1:4).eq.'HG13') then + catm(2,i)(1:4)='2HG1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG13 to 2HG1' + endif + endif + if(catm(2,i)(1:4).eq.' CD1') then + catm(2,i)(1:4)=' CD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'CD1 to CD' + endif + endif + if(catm(2,i)(1:4).eq.' HD1') then + catm(2,i)(1:4)='2HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD1 to 2HD' + endif + endif + if(catm(2,i)(1:4).eq.' HD2') then + catm(2,i)(1:4)='3HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD2 to 3HD' + endif + endif + if(catm(2,i)(1:4).eq.' HD3') then + catm(2,i)(1:4)='4HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD3 to 4HD' + endif + endif + if(catm(2,i)(1:4).eq.'HD11') then + catm(2,i)(1:4)='2HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD11 to 2HD' + endif + endif + if(catm(2,i)(1:4).eq.'HD12') then + catm(2,i)(1:4)='3HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD12 to 3HD' + endif + endif + if(catm(2,i)(1:4).eq.'HD13') then + catm(2,i)(1:4)='4HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD13 to 4HD' + endif + endif + endif + 12 continue +c +c LEU HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c HD11 -> 2HD1 +c HD12 -> 3HD1 +c HD13 -> 4HD1 +c HD21 -> 2HD2 +c HD22 -> 3HD2 +c HD23 -> 4HD2 +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. +c + cres='LEU' + do 13 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'LEU') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.'HD11') then + catm(2,i)(1:4)='2HD1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD11 to 2HD1' + endif + endif + if(catm(2,i)(1:4).eq.'HD12') then + catm(2,i)(1:4)='3HD1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD12 to 3HD1' + endif + endif + if(catm(2,i)(1:4).eq.'HD13') then + catm(2,i)(1:4)='4HD1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD13 to 4HD1' + endif + endif + if(catm(2,i)(1:4).eq.'HD21') then + catm(2,i)(1:4)='2HD2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD21 to 2HD2' + endif + endif + if(catm(2,i)(1:4).eq.'HD22') then + catm(2,i)(1:4)='3HD2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD22 to 3HD2' + endif + endif + if(catm(2,i)(1:4).eq.'HD23') then + catm(2,i)(1:4)='4HD2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD23 to 4HD2' + endif + endif + endif + 13 continue +c +c LYS HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c HG1 -> 2HG +c HG2 -> 3HG +c HG3 -> 2HG +c HD1 -> 2HD +c HD2 -> 3HD +c HD3 -> 2HD +c HE1 -> 2HE +c HE2 -> 3HE +c HE3 -> 2HE +c HZ1 -> 2HZ +c HZ2 -> 3HZ +c HZ3 -> 4HZ +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. Similarly for HG, HD, and HE atoms. +c + cres='LYS' + do 14 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'LYS') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HG1') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG1 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG2') then + catm(2,i)(1:4)='3HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG2 to 3HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG3') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG3 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.' HD1') then + catm(2,i)(1:4)='2HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD1 to 2HD' + endif + endif + if(catm(2,i)(1:4).eq.' HD2') then + catm(2,i)(1:4)='3HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD2 to 3HD' + endif + endif + if(catm(2,i)(1:4).eq.' HD3') then + catm(2,i)(1:4)='2HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD3 to 2HD' + endif + endif + if(catm(2,i)(1:4).eq.' HE1') then + catm(2,i)(1:4)='2HE ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HE1 to 2HE' + endif + endif + if(catm(2,i)(1:4).eq.' HE2') then + catm(2,i)(1:4)='3HE ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HE2 to 3HE' + endif + endif + if(catm(2,i)(1:4).eq.' HE3') then + catm(2,i)(1:4)='2HE ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HE3 to 2HE' + endif + endif + if(catm(2,i)(1:4).eq.' HZ1') then + catm(2,i)(1:4)='2HZ ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HZ1 to 2HZ' + endif + endif + if(catm(2,i)(1:4).eq.' HZ2') then + catm(2,i)(1:4)='3HZ ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HZ2 to 3HZ' + endif + endif + if(catm(2,i)(1:4).eq.' HZ3') then + catm(2,i)(1:4)='4HZ ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HZ3 to 4HZ' + endif + endif + endif + 14 continue +c +c MET HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c HG1 -> 2HG +c HG2 -> 3HG +c HG3 -> 2HG +c HE1 -> 2HE +c HE2 -> 3HE +c HE3 -> 4HE +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. Similarly HG1, HG2. +c + cres='MET' + do 15 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'MET') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HG1') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG1 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG2') then + catm(2,i)(1:4)='3HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG2 to 3HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG3') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG3 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.' HE1') then + catm(2,i)(1:4)='2HE ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HE1 to 2HE' + endif + endif + if(catm(2,i)(1:4).eq.' HE2') then + catm(2,i)(1:4)='3HE ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HE2 to 3HE' + endif + endif + if(catm(2,i)(1:4).eq.' HE3') then + catm(2,i)(1:4)='4HE ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HE3 to 4HE' + endif + endif + endif + 15 continue +c +c PHE HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. +c + cres='PHE' + do 16 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'PHE') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' H1 ') then + catm(2,i)(1:4)='2H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' H2 ') then + catm(2,i)(1:4)='3H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' H3 ') then + catm(2,i)(1:4)='4H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + endif + 16 continue +c +c PRO HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c HG1 -> 2HG +c HG2 -> 3HG +c HG3 -> 2HG +c HD1 -> 2HD +c HD2 -> 3HD +c HD3 -> 2HD +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. +c + cres='PRO' + do 17 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'PRO') then + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HG1') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG1 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG2') then + catm(2,i)(1:4)='3HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG2 to 3HG' + endif + endif + if(catm(2,i)(1:4).eq.' HG3') then + catm(2,i)(1:4)='2HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG3 to 2HG' + endif + endif + if(catm(2,i)(1:4).eq.' HD1') then + catm(2,i)(1:4)='2HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD1 to 2HD' + endif + endif + if(catm(2,i)(1:4).eq.' HD2') then + catm(2,i)(1:4)='3HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD2 to 3HD' + endif + endif + if(catm(2,i)(1:4).eq.' HD3') then + catm(2,i)(1:4)='2HD ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HD3 to 2HD' + endif + endif + endif + 17 continue +c +c SER HN -> H +c H1 -> 2H +c H2 -> 3H +c H3 -> 4H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c HG1 -> HG +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. +c + cres='SER' + do 18 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'SER') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' H1 ') then + catm(2,i)(1:4)='2H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'H1 to 2H' + endif + endif + if(catm(2,i)(1:4).eq.' H2 ') then + catm(2,i)(1:4)='3H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'H2 to 3H' + endif + endif + if(catm(2,i)(1:4).eq.' H3 ') then + catm(2,i)(1:4)='4H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'H3 to 4H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HG1') then + catm(2,i)(1:4)=' HG ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG1 to HG' + endif + endif + endif + 18 continue +c +c THR HN -> H +c HG21 -> 2HG2 +c HG22 -> 3HG2 +c HG23 -> 4HG2 +c + cres='THR' + do 19 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'THR') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.'HG21') then + catm(2,i)(1:4)='2HG2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG21 to 2HG2' + endif + endif + if(catm(2,i)(1:4).eq.'HG22') then + catm(2,i)(1:4)='3HG2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG22 to 3HG2' + endif + endif + if(catm(2,i)(1:4).eq.'HG23') then + catm(2,i)(1:4)='4HG2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG22 to 4HG2' + endif + endif + endif + 19 continue +c +c TRP HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. +c + cres='TRP' + do 20 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'TRP') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + endif + 20 continue +c +c TYR HN -> H +c HB1 -> 2HB +c HB2 -> 3HB +c HB3 -> 2HB +c +c Some PDB files list HB1, HB2, some list HB2, HB3, ultimately both +c cases need to map to 2HB, 3HB. Hence the funky mapping of HB1, +c HB2, and HB3. +c + cres='TYR' + do 21 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'TYR') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' HB1') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB1 to 2HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB2') then + catm(2,i)(1:4)='3HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB2 to 3HB' + endif + endif + if(catm(2,i)(1:4).eq.' HB3') then + catm(2,i)(1:4)='2HB ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HB3 to 2HB' + endif + endif + endif + 21 continue +c +c VAL HN -> H +c H1 -> 2H +c H2 -> 3H +c H3 -> 4H +c HG11 -> 2HG1 +c HG12 -> 3HG1 +c HG13 -> 4HG1 +c HG21 -> 2HG2 +c HG22 -> 3HG2 +c HG23 -> 4HG2 +c + cres='VAL' + do 22 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'VAL') then + if(catm(2,i)(1:4).eq.' HN ') then + catm(2,i)(1:4)=' H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HN to H' + endif + endif + if(catm(2,i)(1:4).eq.' H1 ') then + catm(2,i)(1:4)='2H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'H1 to 2H' + endif + endif + if(catm(2,i)(1:4).eq.' H2 ') then + catm(2,i)(1:4)='3H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'H2 to 3H' + endif + endif + if(catm(2,i)(1:4).eq.' H3 ') then + catm(2,i)(1:4)='4H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'H3 to 4H' + endif + endif + if(catm(2,i)(1:4).eq.'HG11') then + catm(2,i)(1:4)='2HG1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG11 to 2HG1' + endif + endif + if(catm(2,i)(1:4).eq.'HG12') then + catm(2,i)(1:4)='3HG1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG12 to 3HG1' + endif + endif + if(catm(2,i)(1:4).eq.'HG13') then + catm(2,i)(1:4)='4HG1' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG13 to 4HG1' + endif + endif + if(catm(2,i)(1:4).eq.'HG21') then + catm(2,i)(1:4)='2HG2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG21 to 2HG2' + endif + endif + if(catm(2,i)(1:4).eq.'HG22') then + catm(2,i)(1:4)='3HG2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG22 to 3HG2' + endif + endif + if(catm(2,i)(1:4).eq.'HG23') then + catm(2,i)(1:4)='4HG2' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HG22 to 4HG2' + endif + endif + endif + 22 continue +c +c DA O3' -> O3* +c O5' -> O5* +c C5' -> C5* +c H5' -> 2H5* +c H5'' -> 3H5* +c C4' -> C4* +c H4' -> H4* +c O4' -> O4* +c C1' -> C1* +c H1' -> H1* +c C3' -> C3* +c H3' -> H3* +c C2' -> C2* +c H2' -> 2H2* +c H2'' -> 3H2* +c OP1 -> O1P +c OP2 -> O2P +c H61 -> 2H6 +c H62 -> 3H6 +c HO3' -> H3T +c HO5' -> H5T +c + cres='DA' + do 23 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'DA '.or. + + cseq(1,latm(1,i))(1:3).eq.'DA3'.or. + + cseq(1,latm(1,i))(1:3).eq.'DA5') then + if(catm(2,i)(1:4).eq." O3'") then + catm(2,i)(1:4)=' O3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O3' to O3*" + endif + endif + if(catm(2,i)(1:4).eq." O5'") then + catm(2,i)(1:4)=' O5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O5' to O5*" + endif + endif + if(catm(2,i)(1:4).eq." C5'") then + catm(2,i)(1:4)=' C5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C5' to C5*" + endif + endif + if(catm(2,i)(1:4).eq." H5'") then + catm(2,i)(1:4)='2H5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H5' to 2H5*" + endif + endif + if(catm(2,i)(1:4).eq."H5''") then + catm(2,i)(1:4)='3H5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H5'' to 3H5*" + endif + endif + if(catm(2,i)(1:4).eq." C4'") then + catm(2,i)(1:4)=' C4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C4' to C4*" + endif + endif + if(catm(2,i)(1:4).eq." H4'") then + catm(2,i)(1:4)=' H4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H4' to H4*" + endif + endif + if(catm(2,i)(1:4).eq." O4'") then + catm(2,i)(1:4)=' O4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O4' to O4*" + endif + endif + if(catm(2,i)(1:4).eq." C1'") then + catm(2,i)(1:4)=' C1*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C1' to C1*" + endif + endif + if(catm(2,i)(1:4).eq." H1'") then + catm(2,i)(1:4)=' H1*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H1' to H1*" + endif + endif + if(catm(2,i)(1:4).eq." C3'") then + catm(2,i)(1:4)=' C3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C3' to C3*" + endif + endif + if(catm(2,i)(1:4).eq." H3'") then + catm(2,i)(1:4)=' H3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H3' to H3*" + endif + endif + if(catm(2,i)(1:4).eq." C2'") then + catm(2,i)(1:4)=' C2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C2' to C2*" + endif + endif + if(catm(2,i)(1:4).eq." H2'") then + catm(2,i)(1:4)='2H2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H2' to 2H2*" + endif + endif + if(catm(2,i)(1:4).eq."H2''") then + catm(2,i)(1:4)='3H2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H2'' to 3H2*" + endif + endif + if(catm(2,i)(1:4).eq." OP1") then + catm(2,i)(1:4)=' O1P' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"OP1 to O1P" + endif + endif + if(catm(2,i)(1:4).eq." OP2") then + catm(2,i)(1:4)=' O2P' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"OP2 to O2P" + endif + endif + if(catm(2,i)(1:4).eq." H61") then + catm(2,i)(1:4)='2H6 ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H61 to 2H6" + endif + endif + if(catm(2,i)(1:4).eq." H62") then + catm(2,i)(1:4)='3H6 ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H61 to 2H6" + endif + endif + if(catm(2,i)(1:4).eq."HO3'") then + catm(2,i)(1:4)=' H3T' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"HO3' to H3T" + endif + endif + if(catm(2,i)(1:4).eq."HO5'") then + catm(2,i)(1:4)=' H5T' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"HO5' to H5T" + endif + endif + endif + 23 continue +c +c DT O3' -> O3* +c O5' -> O5* +c C5' -> C5* +c H5' -> 2H5* +c H5'' -> 3H5* +c C7 -> C5M +c H71 -> 2H5M +c H72 -> 3H5M +c H73 -> 4H5M +c C4' -> C4* +c H4' -> H4* +c O4' -> O4* +c C1' -> C1* +c H1' -> H1* +c C3' -> C3* +c H3' -> H3* +c C2' -> C2* +c H2' -> 2H2* +c H2'' -> 3H2* +c OP1 -> O1P +c OP2 -> O2P +c HO3' -> H3T +c HO5' -> H5T +c + cres='DT' + do 24 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'DT '.or. + + cseq(1,latm(1,i))(1:3).eq.'DT3'.or. + + cseq(1,latm(1,i))(1:3).eq.'DT5') then + if(catm(2,i)(1:4).eq." O3'") then + catm(2,i)(1:4)=' O3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O3' to O3*" + endif + endif + if(catm(2,i)(1:4).eq." O5'") then + catm(2,i)(1:4)=' O5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O5' to O5*" + endif + endif + if(catm(2,i)(1:4).eq." C5'") then + catm(2,i)(1:4)=' C5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C5' to C5*" + endif + endif + if(catm(2,i)(1:4).eq." H5'") then + catm(2,i)(1:4)='2H5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H5' to 2H5*" + endif + endif + if(catm(2,i)(1:4).eq."H5''") then + catm(2,i)(1:4)='3H5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H5'' to 3H5*" + endif + endif + if(catm(2,i)(1:4).eq." C7 ") then + catm(2,i)(1:4)=' C5M' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C7 to C5M" + endif + endif + if(catm(2,i)(1:4).eq." H71") then + catm(2,i)(1:4)='2H5M' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H71 to 2C5M" + endif + endif + if(catm(2,i)(1:4).eq." H72") then + catm(2,i)(1:4)='3H5M' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H72 to 3C5M" + endif + endif + if(catm(2,i)(1:4).eq." H73") then + catm(2,i)(1:4)='4H5M' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H73 to 4C5M" + endif + endif + if(catm(2,i)(1:4).eq." C4'") then + catm(2,i)(1:4)=' C4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C4' to C4*" + endif + endif + if(catm(2,i)(1:4).eq." H4'") then + catm(2,i)(1:4)=' H4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H4' to H4*" + endif + endif + if(catm(2,i)(1:4).eq." O4'") then + catm(2,i)(1:4)=' O4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O4' to O4*" + endif + endif + if(catm(2,i)(1:4).eq." C1'") then + catm(2,i)(1:4)=' C1*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C1' to C1*" + endif + endif + if(catm(2,i)(1:4).eq." H1'") then + catm(2,i)(1:4)=' H1*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H1' to H1*" + endif + endif + if(catm(2,i)(1:4).eq." C3'") then + catm(2,i)(1:4)=' C3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C3' to C3*" + endif + endif + if(catm(2,i)(1:4).eq." H3'") then + catm(2,i)(1:4)=' H3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H3' to H3*" + endif + endif + if(catm(2,i)(1:4).eq." C2'") then + catm(2,i)(1:4)=' C2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C2' to C2*" + endif + endif + if(catm(2,i)(1:4).eq." H2'") then + catm(2,i)(1:4)='2H2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H2' to 2H2*" + endif + endif + if(catm(2,i)(1:4).eq."H2''") then + catm(2,i)(1:4)='3H2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H2'' to 3H2*" + endif + endif + if(catm(2,i)(1:4).eq." OP1") then + catm(2,i)(1:4)=' O1P' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"OP1 to O1P" + endif + endif + if(catm(2,i)(1:4).eq." OP2") then + catm(2,i)(1:4)=' O2P' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"OP2 to O2P" + endif + endif + if(catm(2,i)(1:4).eq."HO3'") then + catm(2,i)(1:4)=' H3T' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"HO3' to H3T" + endif + endif + if(catm(2,i)(1:4).eq."HO5'") then + catm(2,i)(1:4)=' H5T' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"HO5' to H5T" + endif + endif + endif + 24 continue +c +c DC O3' -> O3* +c O5' -> O5* +c C5' -> C5* +c H5' -> 2H5* +c H5'' -> 3H5* +c C4' -> C4* +c H4' -> H4* +c O4' -> O4* +c C1' -> C1* +c H1' -> H1* +c C3' -> C3* +c H3' -> H3* +c C2' -> C2* +c H2' -> 2H2* +c H2'' -> 3H2* +c HO3' -> H3T +c HO5' -> H5T +c OP1 -> O1P +c OP2 -> O2P +c H41 -> 2H4 +c H42 -> 3H4 +c + cres='DC' + do 25 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'DC '.or. + + cseq(1,latm(1,i))(1:3).eq.'DC3'.or. + + cseq(1,latm(1,i))(1:3).eq.'DC5') then + if(catm(2,i)(1:4).eq." O3'") then + catm(2,i)(1:4)=' O3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O3' to O3*" + endif + endif + if(catm(2,i)(1:4).eq." O5'") then + catm(2,i)(1:4)=' O5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O5' to O5*" + endif + endif + if(catm(2,i)(1:4).eq." C5'") then + catm(2,i)(1:4)=' C5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C5' to C5*" + endif + endif + if(catm(2,i)(1:4).eq." H5'") then + catm(2,i)(1:4)='2H5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H5' to 2H5*" + endif + endif + if(catm(2,i)(1:4).eq."H5''") then + catm(2,i)(1:4)='3H5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H5'' to 3H5*" + endif + endif + if(catm(2,i)(1:4).eq." C4'") then + catm(2,i)(1:4)=' C4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C4' to C4*" + endif + endif + if(catm(2,i)(1:4).eq." H4'") then + catm(2,i)(1:4)=' H4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H4' to H4*" + endif + endif + if(catm(2,i)(1:4).eq." O4'") then + catm(2,i)(1:4)=' O4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O4' to O4*" + endif + endif + if(catm(2,i)(1:4).eq." C1'") then + catm(2,i)(1:4)=' C1*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C1' to C1*" + endif + endif + if(catm(2,i)(1:4).eq." H1'") then + catm(2,i)(1:4)=' H1*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H1' to H1*" + endif + endif + if(catm(2,i)(1:4).eq." C3'") then + catm(2,i)(1:4)=' C3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C3' to C3*" + endif + endif + if(catm(2,i)(1:4).eq." H3'") then + catm(2,i)(1:4)=' H3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H3' to H3*" + endif + endif + if(catm(2,i)(1:4).eq." C2'") then + catm(2,i)(1:4)=' C2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C2' to C2*" + endif + endif + if(catm(2,i)(1:4).eq." H2'") then + catm(2,i)(1:4)='2H2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H2' to 2H2*" + endif + endif + if(catm(2,i)(1:4).eq."H2''") then + catm(2,i)(1:4)='3H2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H2'' to 3H2*" + endif + endif + if(catm(2,i)(1:4).eq."HO3'") then + catm(2,i)(1:4)=' H3T' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"HO3' to H3T" + endif + endif + if(catm(2,i)(1:4).eq."HO5'") then + catm(2,i)(1:4)=' H5T' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"HO5' to H5T" + endif + endif + if(catm(2,i)(1:4).eq." OP1") then + catm(2,i)(1:4)=' O1P' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i)," OP1 to O1P" + endif + endif + if(catm(2,i)(1:4).eq." OP2") then + catm(2,i)(1:4)=' O2P' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i)," OP2 to O2P" + endif + endif + if(catm(2,i)(1:4).eq." H41") then + catm(2,i)(1:4)='2H4 ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i)," H41 to 2H4 " + endif + endif + if(catm(2,i)(1:4).eq." H42") then + catm(2,i)(1:4)='3H4 ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i)," H42 to 3H4 " + endif + endif + endif + 25 continue +c +c DG O3' -> O3* +c O5' -> O5* +c C5' -> C5* +c H5' -> 2H5* +c H5'' -> 3H5* +c C4' -> C4* +c H4' -> H4* +c O4' -> O4* +c C1' -> C1* +c H1' -> H1* +c C3' -> C3* +c H3' -> H3* +c C2' -> C2* +c H2' -> 2H2* +c H2'' -> 3H2* +c OP1 -> O1P +c OP2 -> O2P +c H21 -> 2H2 +c H22 -> 3H2 +c HO3' -> H3T +c HO5' -> H5T +c + cres='DG' + do 26 i=1,natm + if(cseq(1,latm(1,i))(1:3).eq.'DG '.or. + + cseq(1,latm(1,i))(1:3).eq.'DG3'.or. + + cseq(1,latm(1,i))(1:3).eq.'DG5') then + if(catm(2,i)(1:4).eq." O3'") then + catm(2,i)(1:4)=' O3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O3' to O3*" + endif + endif + if(catm(2,i)(1:4).eq." O5'") then + catm(2,i)(1:4)=' O5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O5' to O5*" + endif + endif + if(catm(2,i)(1:4).eq." C5'") then + catm(2,i)(1:4)=' C5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C5' to C5*" + endif + endif + if(catm(2,i)(1:4).eq." H5'") then + catm(2,i)(1:4)='2H5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H5' to 2H5*" + endif + endif + if(catm(2,i)(1:4).eq."H5''") then + catm(2,i)(1:4)='3H5*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H5'' to 3H5*" + endif + endif + if(catm(2,i)(1:4).eq." C4'") then + catm(2,i)(1:4)=' C4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C4' to C4*" + endif + endif + if(catm(2,i)(1:4).eq." H4'") then + catm(2,i)(1:4)=' H4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H4' to H4*" + endif + endif + if(catm(2,i)(1:4).eq." O4'") then + catm(2,i)(1:4)=' O4*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"O4' to O4*" + endif + endif + if(catm(2,i)(1:4).eq." C1'") then + catm(2,i)(1:4)=' C1*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C1' to C1*" + endif + endif + if(catm(2,i)(1:4).eq." H1'") then + catm(2,i)(1:4)=' H1*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H1' to H1*" + endif + endif + if(catm(2,i)(1:4).eq." C3'") then + catm(2,i)(1:4)=' C3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C3' to C3*" + endif + endif + if(catm(2,i)(1:4).eq." H3'") then + catm(2,i)(1:4)=' H3*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H3' to H3*" + endif + endif + if(catm(2,i)(1:4).eq." C2'") then + catm(2,i)(1:4)=' C2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"C2' to C2*" + endif + endif + if(catm(2,i)(1:4).eq." H2'") then + catm(2,i)(1:4)='2H2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H2' to 2H2*" + endif + endif + if(catm(2,i)(1:4).eq."H2''") then + catm(2,i)(1:4)='3H2*' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"H2'' to 3H2*" + endif + endif + if(catm(2,i)(1:4).eq." OP1") then + catm(2,i)(1:4)=' O1P' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i)," OP1 to O1P" + endif + endif + if(catm(2,i)(1:4).eq." OP2") then + catm(2,i)(1:4)=' O2P' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i)," OP2 to O2P" + endif + endif + if(catm(2,i)(1:4).eq." H21") then + catm(2,i)(1:4)='2H2 ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i)," H21 to 2H2" + endif + endif + if(catm(2,i)(1:4).eq." H22") then + catm(2,i)(1:4)='3H2 ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i)," H22 to 3H2" + endif + endif + if(catm(2,i)(1:4).eq."HO3'") then + catm(2,i)(1:4)=' H3T' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"HO3' to H3T" + endif + endif + if(catm(2,i)(1:4).eq."HO5'") then + catm(2,i)(1:4)=' H5T' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),"HO5' to H5T" + endif + endif + endif + 26 continue +c +c Fix N-terminal hydrogens +c +c XAA HT1 -> 2H +c HT2 -> 3H +c HT3 -> 4H +c + do 27 i=1,natm + cres = cseq(1,latm(1,i))(1:3) + if(pre_protein(cseq(1,latm(1,i)))) then + if(catm(2,i)(1:4).eq.' HT1') then + catm(2,i)(1:4)='2H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HT1 to 2H' + endif + endif + if(catm(2,i)(1:4).eq.' HT2') then + catm(2,i)(1:4)='3H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HT2 to 3H' + endif + endif + if(catm(2,i)(1:4).eq.' HT3') then + catm(2,i)(1:4)='4H ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'HT3 to 4H' + endif + endif + endif + 27 continue +c +c Fix C-terminal oxygens +c +c XAA OT1 -> O +c OT2 -> OXT +c + do 28 i=1,natm + cres = cseq(1,latm(1,i))(1:3) + if(pre_protein(cseq(1,latm(1,i)))) then + if(catm(2,i)(1:4).eq.' OT1') then + catm(2,i)(1:4)=' O ' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'OT1 to O' + endif + endif + if(catm(2,i)(1:4).eq.' OT2') then + catm(2,i)(1:4)=' OXT' + if(util_print('sequence',print_high)) then + write(lfnout,2002) cres,latm(1,i),'OT2 to OXT' + endif + endif + endif + 28 continue +c + do 29 i=1,natm + if(latm(2,i).le.0) latm(2,i)=pre_atnum(catm(2,i)(1:2)) + if(latm(2,i).eq.0) then + elemnt(1:1)=' ' + elemnt(2:2)=catm(2,i)(2:2) + latm(2,i)=pre_atnum(elemnt) + if(util_print('sequence',print_high)) then + write(lfnout,2004) latm(1,i) + 2004 format(/,'Unrecognized element in segment ',i4) + endif + endif + 29 continue +c + pre_namiup=.true. +c + return + end diff --git a/contrib/pdb_tools/pdb_amber2nwchem/pre_protein.F b/contrib/pdb_tools/pdb_amber2nwchem/pre_protein.F new file mode 100644 index 0000000000..e411ac8962 --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/pre_protein.F @@ -0,0 +1,35 @@ + logical function pre_protein(sgmnam) +c +c $Id$ +c + implicit none +c +#include "pre_common.fh" +c + character*5 sgmnam +c + pre_protein=.true. +c + if(sgmnam(1:3).eq.'ALA'.or.sgmnam(1:3).eq.'ARG'.or. + + sgmnam(1:3).eq.'ASN'.or.sgmnam(1:3).eq.'ASP'.or. + + sgmnam(1:3).eq.'CYS'.or.sgmnam(1:3).eq.'GLN'.or. + + sgmnam(1:3).eq.'GLU'.or.sgmnam(1:3).eq.'GLY'.or. + + sgmnam(1:3).eq.'HIS'.or.sgmnam(1:3).eq.'ILE'.or. + + sgmnam(1:3).eq.'LEU'.or.sgmnam(1:3).eq.'LYS'.or. + + sgmnam(1:3).eq.'MET'.or.sgmnam(1:3).eq.'PHE'.or. + + sgmnam(1:3).eq.'PRO'.or.sgmnam(1:3).eq.'SER'.or. + + sgmnam(1:3).eq.'THR'.or.sgmnam(1:3).eq.'TRP'.or. + + sgmnam(1:3).eq.'TYR'.or.sgmnam(1:3).eq.'VAL') return + if(ffield(1:5).eq.'amber') then + if(sgmnam(1:3).eq.'ASH'.or.sgmnam(1:3).eq.'CYX'.or. + + sgmnam(1:3).eq.'GLH'.or.sgmnam(1:3).eq.'HID'.or. + + sgmnam(1:3).eq.'HIE'.or.sgmnam(1:3).eq.'HIP') return + endif + if(ffield(1:6).eq.'charmm') then + if(sgmnam(1:3).eq.'CYX'.or.sgmnam(1:3).eq.'HSD'.or. + + sgmnam(1:3).eq.'HSE'.or.sgmnam(1:3).eq.'HSP') return + endif +c + pre_protein=.false. + return + end diff --git a/contrib/pdb_tools/pdb_amber2nwchem/printlevels.fh b/contrib/pdb_tools/pdb_amber2nwchem/printlevels.fh new file mode 100644 index 0000000000..503a77e71f --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/printlevels.fh @@ -0,0 +1,6 @@ + integer print_none, print_low, print_medium, print_high, + $ print_debug, print_default, print_never +C$Id$ + parameter (print_none=0, print_low=10, print_medium=20, + $ print_high=30, print_debug=100, print_never=1000000) + parameter (print_default = print_medium) diff --git a/contrib/pdb_tools/pdb_amber2nwchem/util.fh b/contrib/pdb_tools/pdb_amber2nwchem/util.fh new file mode 100644 index 0000000000..6c01182d3e --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/util.fh @@ -0,0 +1,35 @@ +#include "printlevels.fh" +c +C$Id$ + logical util_print + double precision ddot + double precision util_cpusec + double precision util_wallsec + double precision util_random + integer util_batch_job_time_remaining + integer util_time_remaining + logical util_test_time_remaining + logical util_nwchemrc_get + logical util_module_avail + external util_print + external ddot + external util_cpusec + external util_wallsec + external util_random + external util_batch_job_time_remaining + external util_time_remaining + external util_test_time_remaining + external util_nwchemrc_get + external util_module_avail +c + logical util_xyz_seek + external util_xyz_seek + logical util_xyz_nframes + external util_xyz_nframes +c + logical util_get_io_unit + external util_get_io_unit +c + integer nw_max_path_len ! Maximum path len -> posix standard is what? + parameter (nw_max_path_len = 255) +c diff --git a/contrib/pdb_tools/pdb_amber2nwchem/util_print.F b/contrib/pdb_tools/pdb_amber2nwchem/util_print.F new file mode 100644 index 0000000000..5a1ac017a5 --- /dev/null +++ b/contrib/pdb_tools/pdb_amber2nwchem/util_print.F @@ -0,0 +1,10 @@ + logical function util_print(name, level) + implicit none + character*(*) name + integer level +c + util_print = .false. +c + return +c + end diff --git a/contrib/pdb_tools/pdb_large/GNUmakefile b/contrib/pdb_tools/pdb_large/GNUmakefile new file mode 100644 index 0000000000..c0997872bc --- /dev/null +++ b/contrib/pdb_tools/pdb_large/GNUmakefile @@ -0,0 +1,10 @@ +OBJ_OPTIMIZE=inp.o \ + pdb_large.o \ + pdb_readcommandlinearguments.o + + +all: $(OBJ_OPTIMIZE) + $(FC) -o pdb_large $(OBJ_OPTIMIZE) + +.F.o: + $(FC) -c $< diff --git a/contrib/pdb_tools/pdb_large/README.md b/contrib/pdb_tools/pdb_large/README.md new file mode 100644 index 0000000000..5feafcd316 --- /dev/null +++ b/contrib/pdb_tools/pdb_large/README.md @@ -0,0 +1,12 @@ +# README PDB_LARGE + +NWChem has a special file format for large PDB files. When PDB files are +the residue numbers can no longer be listed with just 4 digits. In the +large format 6 digits can be used. However, for applications it is +useful to be able to convert to and from the large format. + +The large format is specified by providing the PDB keyword "LRGPDB" +before the first "ATOM" or "HETATM" keywork. + +This converter detects what format the input PDB file is in and +generates a PDB output file in the other format. diff --git a/contrib/pdb_tools/pdb_large/inp.F b/contrib/pdb_tools/pdb_large/inp.F new file mode 100644 index 0000000000..b59d52620c --- /dev/null +++ b/contrib/pdb_tools/pdb_large/inp.F @@ -0,0 +1,61 @@ + block data inp_data +C$Id$ + implicit none +#include "inpP.fh" +c + data iread /5/ + data iwrite /6/ + data jrec /-1/ + data jump /0/ + data oswit /.false. / + data nerr /999/ + data nline /0/ + data noline /0/ + data ierrpos /-1/ + data errmsg /' '/ + data input_line /0/ + data xblnk /' '/ +#if defined(SGITFP) || defined(SGI_N32) || defined(CRAY) || defined(HPUX) || defined(WIN32) || defined(PSCALE) || defined(__FLANG) || ( __GNUC__ >= 4) + data xtab /' '/ ! Tab ... no backslash necessary +#elif (defined(LINUX) || defined(MACX)) && !defined(PGLINUX) && !defined(XLFLINUX) &&!( __GNUC__ >= 4) + data xtab /9/ ! Tab ... g77 has trouble with escape sequence +#else + data xtab /'\ '/ ! Tab ... note backslash for cpp +#endif + data xsplit/';'/ + data xcomm /'#'/ + data xback /'\\'/ ! Backslash ... note backslash for cpp + data xquote/'"'/ +c + data save_level /0/ ! Depth of nesting in save/restore + data include_level /0/ +c + end +c +C> \brief Returns the length of the contents of a string +C> +C> In Fortran the length of a string is determined by the amount of +C> memory allocated for it. Often one just wants to know how long the +C> value of the string is without any trailing spaces. Which is what +C> this function returns. +C> +C> \return The length of the value of a string without any trailing +C> spaces +c + integer function inp_strlen(a) + implicit none +#include "inpP.fh" + character*(*) a !< [Input] The string + integer i + logical ois_ws + intrinsic len + character*1 xtest + ois_ws(xtest) = (xtest.eq.xblnk .or. xtest.eq.xtab) +c + do i = len(a),1,-1 + if (.not. ois_ws(a(i:i))) goto 10 + enddo +c + 10 inp_strlen = i +c + end diff --git a/contrib/pdb_tools/pdb_large/inpP.fh b/contrib/pdb_tools/pdb_large/inpP.fh new file mode 100644 index 0000000000..9948bd862a --- /dev/null +++ b/contrib/pdb_tools/pdb_large/inpP.fh @@ -0,0 +1,76 @@ +c +c Private header file for free format input routines +c +** implicit real*8 (a-h,p-w), integer (i-n), logical (o) +** implicit character *8 (z),character *1 (x) +** implicit character *4 (y) +c + integer max_width ! Maximum input width of an input line +C$Id$ + integer max_field ! Maximum no. of fields in an input line + parameter (max_width = 1024) + parameter (max_field = max_width/2 + 1) + character*1024 ja, ia ! Input buffers ... MUST match max_width + character*1024 tmp ! Same size work space + character*80 errmsg ! Error message + character*1 xcomm ! Comment character + character*1 xsplit ! Character to split physical input lines + character*1 xback ! Backslash for concatenation and quoting + character*1 xquote ! Quotation marks for strings + character*1 xblnk ! Space + character*1 xtab ! Tab + common /inp_workc / ia, ja, tmp, errmsg, + $ xcomm, xsplit, xback, xquote, xblnk, xtab +c + integer iread, iwrite + integer jrec ! No. of current field + integer jump ! No. of fields in current line + integer istrt ! Start of fields + integer inumb ! Length of fields + integer nend ! End of fields + integer iwidth ! Length of current logical input line + integer nline ! Current logical line inside physical line + integer noline ! No. of logical lines inside physical line + integer input_line ! No. of current physical input line + integer nerr ! ? + logical oswit ! True if EOF has beeen detected + integer ierrpos ! Input char position where error was detected + integer nstart +c + common /inp_work/ jrec, jump, + $ istrt(max_field), inumb(max_field),iwidth,nend(max_field), + $ nstart(max_field),nline,noline,nerr,oswit,ierrpos, input_line + common/inp_file/iread,iwrite +c + integer max_include_level + parameter (max_include_level = 3) + integer include_level + common/inp_include_info/include_level + character*255 include_file_name + common/inp_include_cinfo/include_file_name(max_include_level) +c +c This drivel below is to save state so that by doing the following +c +c call inp_save_state +c +c call inp_init +c +c call inp_restore_state +c +c It is possible to go off and read another file without +c destroying the state of the current file. Note that this +c only works nested three deep (i.e., 4 files concurrently open) +c +c Note that all of these commons must match in size the +c ones above and that inp_save/restore_state must also be +c consistent +c + character*1 csave + integer isave + integer fsave + integer save_level + common/inp_csave/csave(3*max_width+80+6,3) + common/inp_isave/isave(9+4*max_field,3) + common/inp_fsave/fsave(2,3) + common/inp_save_data/save_level +c diff --git a/contrib/pdb_tools/pdb_large/pdb_large.F b/contrib/pdb_tools/pdb_large/pdb_large.F new file mode 100644 index 0000000000..d383e9c83c --- /dev/null +++ b/contrib/pdb_tools/pdb_large/pdb_large.F @@ -0,0 +1,87 @@ + program pdb_amber2nwchem + implicit none + character(len=1024) :: inputfile + character(len=1024) :: outputfile + character(len=100) :: line + character(len=20) :: l1 ! characters 1-20 + character(len=1) :: chain ! characters 22 + character(len=6) :: resnum ! characters 23-26 + character(len=100) :: l2 ! characters 27-72 + character(len=6) :: old_resnum + integer, external :: inp_strlen + logical, external :: pre_namiup + integer :: lfninp + integer :: lfnout + integer :: iline + integer :: resid + logical :: largepdb + logical :: writtenlarge + call pdb_readcommandlinearguments(inputfile,outputfile) + lfninp = 5 + lfnout = 6 + largepdb = .false. + writtenlarge = .false. + old_resnum = "-1" + resid = 0 + chain = "A" + if (inputfile.ne."-") then + close(lfninp) + open(unit=lfninp,file=inputfile,err=100,form="formatted", + + status="old") + endif + 100 continue + if (outputfile.ne."-") then + close(lfnout) + open(unit=lfnout,file=outputfile,err=200,form="formatted") + endif + 200 continue + iline = 0 + line = "" + do while (line(1:3).ne."END") + iline = iline + 1 + read(lfninp,'(a)')line + select case (line(1:4)) + case ("LRGP") ! LRGPDB + largepdb = .true. + case ("ATOM", "HETA") + if (largepdb) then + read(line,2002,end=300,err=300)l1,resid,l2 + write(lfnout,3002)l1,chain,mod(resid-1,9999)+1, + + l2(1:inp_strlen(l2)) + else + if (.not.writtenlarge) then + write(lfnout,'("LRGPDB")') + writtenlarge = .true. + endif + read(line,1002,end=300,err=300)l1,chain,resnum,l2 + if (resnum.ne.old_resnum) then + resid = resid + 1 + old_resnum = resnum + endif + write(lfnout,2002)l1,resid,l2(1:inp_strlen(l2)) + endif + case ("TER ") + if (largepdb) then + read(line,2002,end=300,err=300)l1,resid,l2 + write(lfnout,3002)l1,chain,mod(resid-1,9999)+1, + + l2(1:inp_strlen(l2)) + else + read(line,1002,end=300,err=300)l1,chain,resnum,l2 + write(lfnout,2002)l1,resid,l2(1:inp_strlen(l2)) + endif + if (chain.eq."Z") then + chain = "A" + else + chain = char(ichar(chain)+1) + endif + case default + write(lfnout,'(a)')line(1:inp_strlen(line)) + end select + enddo + 300 continue + close(lfninp) + close(lfnout) + 1002 format(a20,1x,a1,a4,a) + 2002 format(a20,i6,a) + 3002 format(a20,1x,a1,i4,a) + end diff --git a/contrib/pdb_tools/pdb_large/pdb_readcommandlinearguments.F b/contrib/pdb_tools/pdb_large/pdb_readcommandlinearguments.F new file mode 100644 index 0000000000..f13f1c48f3 --- /dev/null +++ b/contrib/pdb_tools/pdb_large/pdb_readcommandlinearguments.F @@ -0,0 +1,43 @@ + subroutine pdb_readcommandlinearguments(pdb_input,pdb_output) + implicit none + character(len=*), intent(out) :: pdb_input + character(len=*), intent(out) :: pdb_output +! + integer :: narg ! the number of command line arguments + integer :: iarg ! the current command line argument +! + character(len=8) :: arg +! + pdb_input = "-" + pdb_output = "-" + narg = command_argument_count() + select case (narg) + case (2, 4) + iarg=1 + do while (iarg.lt.narg) + call get_command_argument(iarg,arg) + select case (arg) + case ("-i") + iarg=iarg+1 + call get_command_argument(iarg,pdb_input) + case ("-o") + iarg=iarg+1 + call get_command_argument(iarg,pdb_output) + case default + go to 10 + end select + iarg=iarg+1 + enddo + case default + go to 10 + end select + return +! + 10 write(*,*)'Usage: pdb_large -i inputfile -o outputfile' + write(*,*) + write(*,*)'inputfile : PDB file in either normal or large format' + write(*,*) + write(*,*)'outputfile: PDB file in the other format' + stop 10 +! + end diff --git a/contrib/pdb_tools/pdb_scale.F b/contrib/pdb_tools/pdb_scale.F new file mode 100644 index 0000000000..f268e27c92 --- /dev/null +++ b/contrib/pdb_tools/pdb_scale.F @@ -0,0 +1,55 @@ +c======================================================================= +c +!> \brief Scale the coordinates and box size in a PDB file +!> +!> This program scales the atomic coordinates and the box size in a +!> PBD file. The scale factor is provided on the command line. +!> Negative scale factors are allowed but it is not clear what that +!> means if the PDB file has a CRYST1 record. +!> +!> The PDB file is read from standard input and the result written +!> to standard output. +!> + program pdb_scale + character(len=32) :: arg ! command line argument + character(len=80) :: line_in ! line of input file + character(len=80) :: line_out ! line for output file + integer :: iline ! counts lines of input file + integer :: i ! generic counter + double precision :: factor ! scale factor + double precision :: xyz(3) ! coordinates or lattice parameters + iline = 0 + if (command_argument_count().ne.1) then + write(0,*)"ERROR: Invalid invocation of pdb_scale." + write(0,*)" The program expects 1 command line argument" + write(0,*)" that is the scale factor." + endif + call get_command_argument(1,arg) + read(arg,*)factor + do while (.true.) + iline = iline + 1 + read(*,'(a80)',err=999,end=888)line_in + line_out = line_in + select case(line_in(1:6)) + case("ATOM ","HETATM") + read(line_in(31:54),'(3f8.3)')(xyz(i),i=1,3) + do i = 1, 3 + xyz(i) = factor*xyz(i) + enddo + write(line_out(31:54),'(3f8.3)')(xyz(i),i=1,3) + case("CRYST1") + read(line_in(7:33),'(3f9.3)')(xyz(i),i=1,3) + do i = 1, 3 + xyz(i) = factor*xyz(i) + enddo + write(line_out(7:33),'(3f9.3)')(xyz(i),i=1,3) + case default + ! do nothing: the input is echoed to the ouput + end select + write(*,'(a80)')line_out + enddo + 999 write(0,*)"ERROR: I/O error reading line ",iline + 888 continue + end program pdb_scale +c +c======================================================================= diff --git a/contrib/pdb_tools/pdb_supercell.py b/contrib/pdb_tools/pdb_supercell.py new file mode 100644 index 0000000000..717bfdd7b3 --- /dev/null +++ b/contrib/pdb_tools/pdb_supercell.py @@ -0,0 +1,175 @@ +#!/bin/env python +# +# Build a super cell of a given PDB file. +# Because of peculiarities of NWChem we need to produce a +# super cell by creating more atoms with translated coordinates. +# Then we need to write all the non-water atoms out first +# followed by all water atoms. +# +def read_pdb_file(filename): + ''' + Read the contents of the PDB file and return a list of lines. + ''' + fp = open(filename,'r') + pdb = fp.readlines() + fp.close() + return pdb + +def get_cell_dimensions(pdb): + ''' + Take the PDB file as a list of lines. + Read the cell dimensions from the first line (starting with "CRYST1") + and return this as a tuple. + ''' + line = pdb[0] + elements = line.split() + x = float(elements[1]) + y = float(elements[2]) + z = float(elements[3]) + return (x,y,z) + +def gen_translation(direction,pdb): + ''' + Generate the translation vector based on the choice of direction (given as a string), + and the cell dimension of the PDB file. + Return this a tuple. + ''' + (x,y,z) = get_cell_dimensions(pdb) + if direction == "x": + trans = (x,0.0,0.0) + elif direction == "y": + trans = (0.0,y,0.0) + elif direction == "z": + trans = (0.0,0.0,z) + else: + print("Interesting coordinate system you are using. Your choice of direction is %s\n"%direction) + return trans + +def new_cell_dimensions(direction,pdb): + ''' + Given the PDB data and the direction compute the new cell dimensions + after extending the system in the given dimension. + Return the new dimensions as a tuple. + ''' + (x,y,z) = get_cell_dimensions(pdb) + if direction == "x": + cell = (2*x,y,z) + elif direction == "y": + cell = (x,2*y,z) + elif direction == "z": + cell = (x,y,2*z) + else: + print("Interesting coordinate system you are using. Your choice of direction is %s\n"%direction) + return cell + +def split_pdb(pdb): + ''' + Given the PDB file, separate the contents into three parts: + - the cell definition + - the solute atoms + - the solvent atoms + Return these results as a tuple of lists of lines. + ''' + cryst = [] + solute = [] + solvent = [] + for line in pdb: + elements = line.split() + if elements[0] == "CRYST1": + cryst.append(line) + elif elements[0] == "ATOM": + if elements[3] == "HOH": + solvent.append(line) + else: + solute.append(line) + elif elements[0] == "TER": + solute.append(line) + return (cryst,solute,solvent) + +def translate_atoms(atomsin,translation): + ''' + Take a given list of atoms and create a new list of atoms + where all atoms are translated by the given translations. + Return the new list of atoms. The input list remains unchanged. + ''' + (dx,dy,dz) = translation + atomsout = [] + for atom in atomsin: + elements = atom.split() + if elements[0] == "ATOM": + sx = atom[30:38] + sy = atom[38:46] + sz = atom[46:54] + x = float(sx)+dx + y = float(sy)+dy + z = float(sz)+dz + new_atom = "%s%8.3f%8.3f%8.3f%s"%(atom[:30],x,y,z,atom[54:]) + else: + new_atom = atom + atomsout.append(new_atom) + return atomsout + +def new_pdb(cryst,oldsolute,oldsolvent,newsolute,newsolvent): + ''' + Given the different sections of the old PDB and the translated parts of the + solute and solvent lists construct a new PDB file. + Return the list of lines of the new PDB file. + ''' + end = ["END"] + newpdb = cryst+oldsolute+newsolute+oldsolvent+newsolvent+end + return newpdb + +def new_cell(oldcryst,newcell): + ''' + Given the old crystal line and the new cell dimensions construct + a new crystal line. + Return the list of new crystal lines. + ''' + inline = oldcryst[0] + (x,y,z) = newcell + outline = "%s%9.3f%9.3f%9.3f%s"%(inline[:6],x,y,z,inline[33:]) + outlist = [] + outlist.append(outline) + return outlist + +def write_pdb(filename,pdb): + ''' + Given the name of the output file and the list of PDB lines + write the data to the PDB file. + ''' + fp = open(filename,'w') + for line in pdb: + fp.write(line) + fp.close() + +def parse_arguments(): + ''' + Parse command line arguments. + ''' + from argparse import ArgumentParser + prs = ArgumentParser() + prs.add_argument("input",help="the input PDF file") + prs.add_argument("output",help="the output PDF file") + prs.add_argument("direction",help="the direction to expand along") + args = prs.parse_args() + return args + +def execute_with_arguments(args): + inputfile = args.input + outputfile = args.output + direction = args.direction + inputdata = read_pdb_file(inputfile) + trans = gen_translation(direction,inputdata) + newcell = new_cell_dimensions(direction,inputdata) + (cryst,solute,solvent) = split_pdb(inputdata) + newcryst = new_cell(cryst,newcell) + newsolute = translate_atoms(solute,trans) + newsolvent = translate_atoms(solvent,trans) + outputdata = new_pdb(newcryst,solute,solvent,newsolute,newsolvent) + write_pdb(outputfile,outputdata) + +def main(): + execute_with_arguments(parse_arguments()) + +if __name__ == "__main__": + main()