comp-chem/NWChem
1
0
Fork 0
mirror of https://github.com/nwchemgit/nwchem.git synced 2026-07-28 14:15:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

HvD: Adding tools to help with manipulating PDB files.

Browse source
This commit is contained in:
hjjvandam 2020-10-13 13:42:10 -04:00
parent f3ab5e5fb5
commit 98a16a393c
23 changed files with 3177 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

7
contrib/pdb_tools/GNUmakefile Normal file
View file

@ -0,0 +1,7 @@
OBJ_OPTIMIZE=pdb_scale.o
all: $(OBJ_OPTIMIZE)
$(FC) -o pdb_scale $(OBJ_OPTIMIZE)
.F.o:
$(FC) -c $<

33
contrib/pdb_tools/README.md Normal file
View file

@ -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.

14
contrib/pdb_tools/pdb_amber2nwchem/GNUmakefile Normal file
View file

@ -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 $<

38
contrib/pdb_tools/pdb_amber2nwchem/README.md Normal file
View file

@ -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)

61
contrib/pdb_tools/pdb_amber2nwchem/inp.F Normal file
View file

@ -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

76
contrib/pdb_tools/pdb_amber2nwchem/inpP.fh Normal file
View file

@ -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

101
contrib/pdb_tools/pdb_amber2nwchem/pdb_amber2nwchem.F Normal file
View file

@ -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

64
contrib/pdb_tools/pdb_amber2nwchem/pdb_readcommandlinearguments.F Normal file
View file

@ -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

51
contrib/pdb_tools/pdb_amber2nwchem/pre_atnum.F Normal file
View file

@ -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

50
contrib/pdb_tools/pdb_amber2nwchem/pre_common.fh Normal file
View file

@ -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

2077
contrib/pdb_tools/pdb_amber2nwchem/pre_namiup.F Normal file
View file

File diff suppressed because it is too large Load diff

35
contrib/pdb_tools/pdb_amber2nwchem/pre_protein.F Normal file
View file

@ -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

6
contrib/pdb_tools/pdb_amber2nwchem/printlevels.fh Normal file
View file

@ -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)

35
contrib/pdb_tools/pdb_amber2nwchem/util.fh Normal file
View file

@ -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

10
contrib/pdb_tools/pdb_amber2nwchem/util_print.F Normal file
View file

@ -0,0 +1,10 @@
logical function util_print(name, level)
implicit none
character*(*) name
integer level
c
util_print = .false.
c
return
c
end

10
contrib/pdb_tools/pdb_large/GNUmakefile Normal file
View file

@ -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 $<

12
contrib/pdb_tools/pdb_large/README.md Normal file
View file

@ -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.

61
contrib/pdb_tools/pdb_large/inp.F Normal file
View file

@ -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

76
contrib/pdb_tools/pdb_large/inpP.fh Normal file
View file

@ -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

87
contrib/pdb_tools/pdb_large/pdb_large.F Normal file
View file

@ -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

43
contrib/pdb_tools/pdb_large/pdb_readcommandlinearguments.F Normal file
View file

@ -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

55
contrib/pdb_tools/pdb_scale.F Normal file
View file

@ -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=======================================================================

175
contrib/pdb_tools/pdb_supercell.py Normal file
View file

@ -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()