Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F66098492
charmm2lammps.pl
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Jun 8, 06:39
Size
70 KB
Mime Type
text/x-perl
Expires
Mon, Jun 10, 06:39 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18171947
Attached To
rLAMMPS lammps
charmm2lammps.pl
View Options
#!/usr/bin/perl
#
# program: charmm2lammps.pl
# author: Pieter J. in 't Veld,
# pjintve@sandia.gov, veld@verizon.net
# date: February 12-23, April 5, 2005.
# purpose: Translation of charmm input to lammps input
#
# Notes: Copyright by author for Sandia National Laboratories
# 20050212 Needed (in the same directory):
# - $project.crd ; Assumed to be correct and running
# - $project.psf ; CHARMM configs
# - top_$forcefield.rtf ;
# - par_$forcefield.prm ;
# Ouput:
# - $project.data ; LAMMPS data file
# - $project.in ; LAMMPS input file
# - $project_ctrl.pdb ; PDB control file
# - $project_ctrl.psf ; PSF control file
# 20050218 Optimized for memory usage
# 20050221 Rotation added
# 20050222 Water added
# 20050223 Ions added
# 20050405 Water bug fixed; addition of .pdb input
# 20050407 project_ctrl.psf bug fixed; addition of -border
# 20050519 Added interpretation of charmm xplor psfs
# 20050603 Fixed centering issues
# 20050630 Fixed symbol issues arising from salt addition
# 20060818 Changed reading of pdb format to read exact columns
# 20070109 Changed AddMass() to use $max_id correctly
# 20160114 Added compatibility for parameter files that use IMPROPERS instead of IMPROPER
# Print warning when not all parameters are detected. Set correct number of atom types.
# 20160613 Fix off-by-one issue in atom type validation check
# Replace -charmm command line flag with -nohints flag
# and enable type hints in data file by default.
# Add hints also to section headers
# Add a brief minimization to example input template.
# 20161001 Added 'CMAP crossterms' section at the end of the data file
# 20161001 Added instructions in CMAP section to fix problem if 'ter'
# is not designated in the .pdb file to identify last amino acid
# 20161005 Added tweak to embed command line in generated LAMMPS input
#
# General Many thanks to Paul S. Crozier for checking script validity
# against his projects.
# Also thanks to Xiaohu Hu (hux2@ornl.gov) and Robert A. Latour
# (latourr@clemson.edu), David Hyde-Volpe, and Tigran Abramyan,
# Clemson University and Chris Lorenz (chris.lorenz@kcl.ac.uk),
# King's College London for their efforts to add CMAP sections,
# which is implemented using the option flag "-cmap".
# Initialization
sub Test
{
my $name = shift(@_);
printf("Error: file %s not found\n", $name) if (!scalar(stat($name)));
return !scalar(stat($name));
}
sub Initialize # initialization
{
my $k = 0;
my @dir = ("x", "y", "z");
my @options = ("-help", "-nohints", "-water", "-ions", "-center",
"-quiet", "-pdb_ctrl", "-l", "-lx", "-ly", "-lz",
"-border", "-ax", "-ay", "-az", "-cmap");
my @remarks = ("display this message",
"do not print type and style hints in data file",
"add TIP3P water [default: 1 g/cc]",
"add (counter)ions using Na+ and Cl- [default: 0 mol/l]",
"recenter atoms",
"do not print info",
"output project_ctrl.pdb [default: on]",
"set x-, y-, and z-dimensions simultaneously",
"x-dimension of simulation box",
"y-dimension of simulation box",
"z-dimension of simulation box",
"add border to all sides of simulation box [default: 0 A]",
"rotation around x-axis",
"rotation around y-axis",
"rotation around z-axis",
"generate a CMAP section in data file"
);
my $notes;
$program = "charmm2lammps";
$version = "1.9.1";
$year = "2016";
$add = 1;
$water_dens = 0;
$ions = 0;
$info = 1;
$center = 0;
$net_charge = 0;
$ion_molar = 0;
$pdb_ctrl = 1;
$border = 0;
$L = (0, 0, 0);
$cmap = 0;
@R = M_Unit();
$notes = " * The average of extremes is used as the origin\n";
$notes .= " * Residues are numbered sequentially\n";
$notes .= " * Water is added on an FCC lattice: allow 5 ps for";
$notes .= " equilibration\n";
$notes .= " * Ions are added randomly and only when water is present\n";
$notes .= " * CHARMM force field v2.7 parameters used for";
$notes .= " water and NaCl\n";
$notes .= " * Rotation angles are in degrees\n";
$notes .= " * Rotations are executed consecutively: -ax -ay != -ay -ax\n";
$notes .= " * CHARMM files needed in execution directory:\n";
$notes .= " - project.crd coordinates\n";
$notes .= " - project.pdb when project.crd is absent\n";
$notes .= " - project.psf connectivity\n";
$notes .= " - top_forcefield.rtf topology\n";
$notes .= " - par_forcefield.prm parameters\n";
$notes .= " * Output files written to execution directory:\n";
$notes .= " - project.data LAMMPS data file\n";
$notes .= " - project.in suggested LAMMPS input script\n";
$notes .= " - project_ctrl.pdb control file when requested\n";
# record full command line for later use
$cmdline = "$program.pl " . join(" ",@ARGV);
foreach (@ARGV)
{
if (substr($_, 0, 1) eq "-")
{
my $k = 0;
my @tmp = split("=");
my $switch = ($arg[1] eq "")||($arg[1] eq "on")||($arg[1]!=0);
$tmp[0] = lc($tmp[0]);
foreach (@options)
{
last if ($tmp[0] eq substr($_, 0 , length($tmp[0])));
++$k;
}
$help = 1 if (!$k--); # -help
$add = 0 if (!$k--); # -nohints
$water_dens = ($tmp[1] ne "" ? $tmp[1] : 1) if (!$k--); # -water
$ion_molar = abs($tmp[1]) if (!$k); # -ions
$ions = 1 if (!$k--); # ...
$center = 1 if (!$k--); # -center
$info = 0 if (!$k--); # -quiet
$pdb_ctrl = $switch if (!$k--); # -pdb_ctrl
my $flag = $k--; # -l
$L[0] = abs($tmp[1]) if (!($flag && $k--)); # -lx
$L[1] = abs($tmp[1]) if (!($flag && $k--)); # -ly
$L[2] = abs($tmp[1]) if (!($flag && $k--)); # -lz
$border = abs($tmp[1]) if (!$k--); # -border
@R = M_Dot(M_Rotate(0, $tmp[1]), @R) if (!$k--);# -ax
@R = M_Dot(M_Rotate(1, $tmp[1]), @R) if (!$k--);# -ay
@R = M_Dot(M_Rotate(2, $tmp[1]), @R) if (!$k--);# -az
$cmap = ($tmp[1] ne "" ? $tmp[1] : 22) if (!$k--); # -cmap
print("Warning: ignoring unknown command line flag: $tmp[0]\n") unless $k;
}
else
{
$forcefield = $_ if (!$k);
$project = $_ if ($k++ == 1);
}
}
$water_dens = 1 if ($ions && !$water_dens);
if (($k<2)||$help)
{
printf("\n%s v%s (c)2005-%s by Pieter J. in \'t Veld and others\n\n",
$program, $version, $year);
printf("Usage:\n %s.pl [-option[=#] ..] forcefield project\n\n",$program);
printf("Options:\n");
for (my $i=0; $i<scalar(@options); ++$i)
{
printf(" %-10.10s %s\n", $options[$i], $remarks[$i]);
}
printf("\nNotes:\n%s\n", $notes);
exit(-1);
}
else { printf("\n%s v%s (c)2005-%s\n\n", $program, $version, $year) if ($info); }
my $flag = Test($Parameters = "par_$forcefield.prm");
$flag |= Test($Topology = "top_$forcefield.rtf");
$flag |= Test($Pdb = "$project.pdb")
if (!scalar(stat($Crd = "$project.crd")));
$flag |= Test($Psf = "$project.psf") if ($look eq "");
$pdb = ($Pdb ne "") ? 1 : 0;
printf("Conversion aborted\n\n") if ($flag);
exit(-1) if ($flag);
printf("Info: using $Pdb instead of $Crd\n") if (!scalar(stat($Crd)));
for (my $i=0; $i<3; ++$i)
{
printf("Info: l%s not set: will use extremes\n",
("x", "y", "z")[$i]) if ($info&&!$L[$i]);
}
open(PARAMETERS, "<par_$forcefield.prm");
}
# Vector manipulation
sub V_String
{
my @v = @_;
return "{".$v[0].", ".$v[1].", ".$v[2]."}";
}
sub V_Add
{
my @v1 = splice(@_, 0, 3);
my @v2 = splice(@_, 0, 3);
return ($v1[0]+$v2[0], $v1[1]+$v2[1], $v1[2]+$v2[2]);
}
sub V_Subtr
{
my @v1 = splice(@_, 0, 3);
my @v2 = splice(@_, 0, 3);
return ($v1[0]-$v2[0], $v1[1]-$v2[1], $v1[2]-$v2[2]);
}
sub V_Dot
{
my @v1 = splice(@_, 0, 3);
my @v2 = splice(@_, 0, 3);
return $v1[0]*$v2[0]+$v1[1]*$v2[1]+$v1[2]*$v2[2];
}
sub V_Mult
{
my @v = splice(@_, 0, 3);
my $f = shift(@_);
return ($f*$v[0], $f*$v[1], $f*$v[2]);
}
sub M_String
{
my $string;
for (my $i=0; $i<3; ++$i)
{
$string .= ", " if ($i);
$string .= V_String(splice(@_, 0, 3));
}
return "{".$string."}";
}
sub M_Transpose
{
return
(@_[0], @_[3], @_[6],
@_[1], @_[4], @_[7],
@_[2], @_[5], @_[8]);
}
sub M_Dot
{
my @v11 = splice(@_, 0, 3);
my @v12 = splice(@_, 0, 3);
my @v13 = splice(@_, 0, 3);
my @m = M_Transpose(splice(@_, 0, 9));
my @v21 = splice(@m, 0, 3);
my @v22 = splice(@m, 0, 3);
my @v23 = splice(@m, 0, 3);
return (
V_Dot(@v11, @v21), V_Dot(@v11, @v22), V_Dot(@v11, @v23),
V_Dot(@v12, @v21), V_Dot(@v12, @v22), V_Dot(@v12, @v23),
V_Dot(@v13, @v21), V_Dot(@v13, @v22), V_Dot(@v13, @v23));
}
sub M_Unit { return (1,0,0, 0,1,0, 0,0,1); }
sub PI { return 4*atan2(1,1); }
sub M_Rotate
{ # vmd convention
my $n = shift(@_);
my $alpha = shift(@_)*PI()/180;
my $cos = cos($alpha);
my $sin = sin($alpha);
$cos = 0 if (abs($cos)<1e-16);
$sin = 0 if (abs($sin)<1e-16);
return (1,0,0, 0,$cos,-$sin, 0,$sin,$cos) if ($n==0); # around x-axis
return ($cos,0,$sin, 0,1,0, -$sin,0,$cos) if ($n==1); # around y-axis
return ($cos,-$sin,0, $sin,$cos,0, 0,0,1) if ($n==2); # around z-axis
return M_Unit();
}
sub MV_Dot
{
my @v11 = splice(@_, 0, 3);
my @v12 = splice(@_, 0, 3);
my @v13 = splice(@_, 0, 3);
my @v2 = splice(@_, 0, 3);
return (V_Dot(@v11, @v2), V_Dot(@v12, @v2), V_Dot(@v13, @v2));
}
# CHARMM input
sub PSFConnectivity
{
my $n = PSFGoto(bonds);
return if (scalar(@nconnect));
printf("Info: creating connectivity\n") if ($info);
for (my $i=0; $i<$n; ++$i)
{
my @bond = PSFGet(2);
$connect[$bond[0]][$nconnect[$bond[0]]++] = $bond[1];
$connect[$bond[1]][$nconnect[$bond[1]]++] = $bond[0];
}
}
sub PSFDihedrals # hack to accomodate
{ # LAMMPS' way of calc
$idihedral = 0; # LJ 1-4 interactions
return $ndihedral if (($dihedral_flag = $ndihedral ? 1 : 0));
PSFConnectivity();
printf("Info: creating dihedrals\n") if ($info);
my $n = scalar(@nconnect);
my @bonded = ();
for (my $i=1; $i<=$n; ++$i)
{
$bonded[0] = $i;
for (my $i=0; $i<scalar($nconnect[$bonded[0]]); ++$i)
{
$bonded[1] = $connect[$bonded[0]][$i];
for (my $i=0; $i<scalar($nconnect[$bonded[1]]); ++$i)
{
next if (($bonded[2] = $connect[$bonded[1]][$i])==$bonded[0]);
for (my $i=0; $i<scalar($nconnect[$bonded[2]]); ++$i)
{
next if (($bonded[3] = $connect[$bonded[2]][$i])==$bonded[1]);
next if ($bonded[3]<$bonded[0]);
$dihedral[$ndihedral++] = join(" ", @bonded);
}
}
}
}
$dihedral_flag = 1;
return $ndihedral;
}
sub CreatePSFIndex # make an index of id
{ # locations
my @psf_ids = ("!NATOM","!NBOND:","!NTHETA:","!NPHI:","!NIMPHI:");
my @ids = (atoms, bonds, angles, dihedrals, impropers);
my $k = 0;
my %hash;
printf("Info: creating PSF index\n") if ($info);
open(PSF, "<$project.psf") if (fileno(PSF) eq "");
foreach (@psf_ids) { $hash{$_} = shift(@ids); };
while (<PSF>)
{
chop();
my @tmp = split(" ");
my $n = $hash{$tmp[1]};
$PSFIndex{$n} = tell(PSF)." ".$tmp[0] if ($n ne "");
}
}
sub PSFGoto # goto $ident in <PSF>
{
CreatePSFIndex() if (!scalar(%PSFIndex));
my $id = shift(@_);
my @n = split(" ", $PSFIndex{$id});
@PSFBuffer = ();
# return PSFDihedrals() if ($id eq "dihedrals");
if (!scalar(@n))
{
printf("Warning: PSF index for $id not found\n");
seek(PSF, 0, SEEK_END);
return -1;
}
seek(PSF, $n[0], SEEK_SET);
return $n[1];
}
sub PSFGet
{
if ($dihedral_flag)
{
$dihedral_flag = $idihedral+1<$ndihedral ? 1 : 0;
return split(" ", $dihedral[$idihedral++]);
}
if (!scalar(@PSFBuffer))
{
my $line = <PSF>;
chop($line);
@PSFBuffer = split(" ", $line);
}
return splice(@PSFBuffer, 0, shift(@_));
}
sub PSFWrite
{
my $items = shift(@_);
my $n = $items;
if ($psf_ncols>7) { printf(PSF_CTRL "\n"); $psf_ncols = 0; }
foreach(@_)
{
printf(PSF_CTRL " %7d", $_);
++$psf_ncols;
if ((!--$n) && ($psf_ncols>7))
{
printf(PSF_CTRL "\n");
$psf_ncols = 0;
$n = $items;
}
}
}
sub CRDGoto
{
my $n;
return if (shift(@_) ne "atoms");
open(CRD, "<".($pdb ? $Pdb : $Crd)) if (fileno(CRD) eq "");
seek(CRD, 0, SEEK_SET);
return PSFGoto(atoms) if ($pdb);
while (substr($n = <CRD>, 0, 1) eq "*") {}
chop($n);
return $n;
}
sub NextPDB2CRD
{
my @n = (6,5,1,4,1,3,1,1,4,1,3,8,8,8,6,6,6,4,2,2);
my @data = ();
my $c = 0;
my $line;
while (substr($line = <CRD>, 0, 4) ne "ATOM") {};
chop($line);
foreach (@n) { push(@data, substr($line, ($c += $_)-$_, $_)); }
return @data[1, 8, 5, 3, 11, 12, 13, 17, 8, 15];
}
sub Delete
{
my $item = shift(@_);
my $k = 0;
my @list;
foreach (@_)
{
my @tmp = split(" ");
delete($tmp[$item]);
$list[$k++] = join(" ", @tmp);
}
return @list;
}
sub CreateID # create id from list
{
my $n = scalar(@_);
my @list = @_;
my $id = "";
my $flag = $list[0] gt $list[-1];
my $j = $n;
my $tmp;
return "" if (scalar(@list)<$n);
$flag = $list[1] gt $list[-2]
if ((scalar(@list)>3)&&($list[0] eq $list[-1]));
for (my $i=0; $i<$n; ++$i)
{
$id .= ($i ? " " : "").($tmp = $list[$flag ? --$j : $i]);
$id .= substr(" ", 0, 4-length($tmp));
}
return $id;
}
sub AtomTypes
{
my $n = PSFGoto(atoms);
my %list;
return () if ($n<1);
$atom_types[0] = -1;
for (my $i=0; $i<$n; ++$i)
{
my @tmp = split(" ", <PSF>);
$tmp[5] = $symbols{$tmp[5]}
if ((substr($tmp[5],0,1) lt '0')||(substr($tmp[5],0,1) gt '9'));
push(@atom_types, $tmp[5]);
++$list{$tmp[5]};
}
if ($water_dens)
{
push(@atom_types, $symbols{HT}); ++$list{$symbols{HT}};
push(@atom_types, $symbols{OT}); ++$list{$symbols{OT}};
}
if ($ions)
{
push(@atom_types, $symbols{CLA}); ++$list{$symbols{CLA}};
push(@atom_types, $symbols{SOD}); ++$list{$symbols{SOD}};
}
return sort({$a<=>$b} keys(%list));
}
sub Markers
{
my %markers = (
NONBONDED => '0',
BONDS => '1',
ANGLES => '2',
DIHEDRALS => '3',
IMPROPERS => '4',
IMPROPER => '4'
);
return %markers;
}
sub NonBond
{
my @cols = @_;
my $f = (scalar(@cols)>3)&&(substr($cols[3],0,1) ne "!");
my @tmp = (-$cols[1], $cols[2],
$f ? -$cols[4]:-$cols[1], $f ? $cols[5]:$cols[2]);
$tmp[1] *= 2.0**(5/6); # adjust sigma
$tmp[3] *= 2.0**(5/6); # adjust sigma 1-4
return join(" ", @tmp);
}
sub AtomParameters # non-bonded parameters
{
my @types;
my @list;
my $k = 0;
my $read = 0;
my %markers = Markers();
foreach(@_) { $types{$ids{$_}} = $k++; }
seek(PARAMETERS, 0, 0);
while (<PARAMETERS>)
{
chop();
my @cols = split(" ");
if ($read&&(scalar(@cols)>1)&&
(substr($cols[0],0,1) ne "!")&&($cols[1] lt "A"))
{
my $k = $types{shift(@cols)};
$list[$k] = NonBond(@cols) if ($k ne "");
}
if ($markers{$cols[0]} ne "") {
$read = ($markers{$cols[0]} eq "0") ? 1 : 0; }
}
$list[$types{HT}] = NonBond(0, -0.046, 0.2245)
if ($water_dens&&($list[$types{HT}] eq ""));
$list[$types{OT}] = NonBond(0, -0.152100, 1.768200)
if ($water_dens&&($list[$types{OT}] eq ""));
$list[$types{CLA}] = NonBond(0, -0.150, 2.27)
if ($ions&&($list[$types{CLA}] eq ""));
$list[$types{SOD}] = NonBond(0, -0.0469, 1.36375)
if ($ions&&($list[$types{SOD}] eq ""));
return @list;
}
sub BondedTypes # create bonded types
{
my $mode = shift(@_); # operation mode
my $items = (2, 3, 4, 4)[$mode]; # items per entry
my $id = (bonds, angles, dihedrals, impropers)[$mode];
my $n = PSFGoto($id);
my %list;
for (my $i=0; $i<$n; ++$i)
{
my @tmp = ();
foreach (PSFGet($items)) { push(@tmp, $ids{$atom_types[$_]}); }
++$list{CreateID(@tmp)};
}
++$list{CreateID(HT, OT)} if ($water_dens&&($mode==0));
++$list{CreateID(HT, OT, HT)} if ($water_dens&&($mode==1));
@types = sort(keys(%list));
}
sub Parameters # parms from columns
{
my $items = shift(@_);
my @cols = @_;
my $parms = "";
for (my $i=$items; ($i<scalar(@cols))&&(substr($cols[$i],0,1)ne"!"); ++$i)
{
$parms = $parms.($i>$items ? " " : "").$cols[$i];
}
return $parms;
}
sub BondedParameters # distil parms from
{ # <PARAMETERS>
my $mode = shift(@_); # bonded mode
return if (($mode>3)||($mode<0));
my $items = (2, 3, 4, 4)[$mode]; # items per entry
my $name = ("bond", "angle", "dihedral", "improper")[$mode];
my $read = 0;
my $k = 0;
my %markers = Markers();
my @set;
my @tmp;
my $f;
my %list;
my %link;
@parms = ();
foreach(@types) { $link{$_} = $k++; }
seek(PARAMETERS, 0, 0);
while (<PARAMETERS>)
{
chomp();
my @cols = split(" ");
if ($read&&(scalar(@cols)>$items)&&($cols[$items] lt "A"))
{
if (($items==4)&&(($f = ($cols[1] eq "X")&&($cols[2] eq "X"))||
(($cols[0] eq "X")&&($cols[3] eq "X")))) # wildcards
{
my $id = CreateID(($cols[1-$f], $cols[2+$f]));
for ($k=0; $k<scalar(@types); ++$k)
{
if (!$set[$k])
{
my @tmp = split(" ", $types[$k]);
if (CreateID($tmp[1-$f], $tmp[2+$f]) eq $id)
{
if ($mode==2)
{
if ($parms[$k] eq "") {
$parms[$k] = Parameters($items,@cols)." 1"; }
else {
$parms[$k] .= ":".Parameters($items,@cols)." 0"; }
}
else {
$parms[$k] .= Parameters($items,@cols); }
}
}
}
}
else # regular
{
for (my $i=0; $i<$items; ++$i) { $tmp[$i] = $cols[$i]; };
$k = $link{CreateID(@tmp)};
if ($k ne "")
{
$parms[$k] = "" if (!$set[$k]);
$parms[$k] .= ($set[$k]++ ? ":" : "").Parameters($items,@cols);
$parms[$k] .= ($set[$k]-1 ? " 0" : " 1") if ($mode==2);
}
}
}
if ($markers{$cols[0]}) {
$read = ($markers{$cols[0]} eq $mode+1) ? 1 : 0; }
}
if ($water_dens)
{
$parms[$link{CreateID(HT, OT)}] = "450 0.9572" if ($mode==0);
$parms[$link{CreateID(HT, OT, HT)}] = "55 104.52" if ($mode==1);
}
for (my $i=0; $i<scalar(@types); ++$i)
{
printf("Warning: %s parameter %4d for [%s] was not found\n",
$name, $i+1, $types[$i]) if ($parms[$i] eq "");
}
}
sub SetScreeningFactor # set screening factor
{
my $id = shift(@_);
my $value = shift(@_);
my $new = "";
foreach (split(":", $parms[$id]))
{
my @tmp = split(" ");
$tmp[-1] = $value if ($tmp[-1]);
$new .= ":" if ($new ne "");
$new .= join(" ", @tmp);
}
$parms[$id] = $new;
}
sub CorrectDihedralParameters
{
my $n = PSFGoto(dihedrals);
my %hash;
my $hash_id;
my $id1;
my $id2;
my $first;
my $last;
for (my $i=0; $i<$n; ++$i)
{
my @bonded = PSFGet(4);
my @tmp = ();
foreach (@bonded) { push(@tmp, $ids{$atom_types[$_]}); }
$id1 = $link{CreateID(@tmp)}-1;
$first = $bonded[0];
$last = $bonded[3];
if ($first>$last) { my $tmp = $first; $first = $last; $last = $tmp; }
if (($id2 = $hash{$hash_id = $first." ".$last}) eq "")
{
$hash{$hash_id} = $id1; # add id to hash
}
else
{
SetScreeningFactor($id1, 0.5); # 6-ring: shared 1-4
SetScreeningFactor($id2, 0.5);
}
}
$n = PSFGoto(angles);
for (my $i=0; $i<$n; ++$i)
{
my @bonded = PSFGet(3);
$first = $bonded[0];
$last = $bonded[2];
if ($first>$last) { my $tmp = $first; $first = $last; $last = $tmp; }
if (($id1 = $hash{$first." ".$last}) ne "")
{
SetScreeningFactor($id1, 0); # 5-ring: no 1-4
}
}
}
sub AddMass
{
my $symbol = shift(@_);
my $mass = shift(@_);
return if ($symbols{$symbol} ne "");
$ids{++$max_id} = $symbol;
$masses{$max_id} = $mass;
$symbols{$symbol} = $max_id;
}
sub ReadTopology # read topology links
{
my $id = shift(@_);
my $item = shift(@_);
my $read = 0;
my @tmp;
open(TOPOLOGY, "<top_$forcefield.rtf");
$max_id = 0;
while (<TOPOLOGY>)
{
chop(); # delete CR at end
my @tmp = split(" ");
$read = 1 if ($tmp[0] eq "MASS");
if ($read&&($tmp[0] eq "MASS"))
{
$symbols{$tmp[2]} = $tmp[1];
$ids{$tmp[1]} = $tmp[2];
$masses{$tmp[1]} = $tmp[3];
$max_id = $tmp[1] if ($max_id<$tmp[1]);
}
# $names{$tmp[1]} = $tmp[4] if ($read&&($tmp[0] eq "MASS"));
last if ($read&&!scalar(@tmp)); # quit reading
}
AddMass(HT, 1.00800);
AddMass(OT, 15.99940);
AddMass(CLA, 35.450000);
AddMass(SOD, 22.989770);
close(TOPOLOGY);
}
sub CrossLink # symbolic cross-links
{
my @list = @_;
my $n = scalar(@list);
my %hash;
for (my $i=0; $i<$n; ++$i) { $hash{$list[$i]} = $i+1; }
return %hash;
}
sub CharacterizeBox
{
my $flag = 1;
my @x = (-$L[0]/2, $L[0]/2);
my @y = (-$L[1]/2, $L[1]/2);
my @z = (-$L[2]/2, $L[2]/2);
my $n = CRDGoto(atoms);
my $extremes = !($L[0] && $L[1] && $L[2]);
@Center = (0, 0, 0);
return if (!$n);
for (my $i=0; $i<$n; ++$i)
{
my @tmp = $pdb ? NextPDB2CRD() : split(" ", <CRD>);
my @p = @tmp[-6, -5, -4];
@p = MV_Dot(@R, @p);
$x[0] = $p[0] if ($flag||($p[0]<$x[0]));
$x[1] = $p[0] if ($flag||($p[0]>$x[1]));
$y[0] = $p[1] if ($flag||($p[1]<$y[0]));
$y[1] = $p[1] if ($flag||($p[1]>$y[1]));
$z[0] = $p[2] if ($flag||($p[2]<$z[0]));
$z[1] = $p[2] if ($flag||($p[2]>$z[1]));
$flag = 0 if ($flag);
}
$L[0] = $x[1]-$x[0] if (!$L[0]);
$L[1] = $y[1]-$y[0] if (!$L[1]);
$L[2] = $z[1]-$z[0] if (!$L[2]);
$L[0] += $border;
$L[1] += $border;
$L[2] += $border;
@Center = (($x[1]+$x[0])/2, ($y[1]+$y[0])/2, ($z[1]+$z[0])/2);
printf("Info: recentering atoms\n") if ($info&&$center);
}
sub SetupWater
{
return if (!$water_dens);
my $dens = 1000*$water_dens; # kg/m^3
my $m = 0.018; # kg/mol
my $loh = 0.9572; # l[O-H] in [A]
my $s_OT = 1.7682; # CHARMM sigma [A]
my $ahoh = (180-104.52)/360*PI();
my @p = ($loh*cos($ahoh), $loh*sin($ahoh), 0);
printf("Info: creating fcc water\n") if ($info);
$n_water = 4; # molecules/cell
$nav = 6.022e23; # 1/mol
$v_water = $m/$nav/$dens*1e30; # water volume [A^3]
$r_water = $s_OT*2**(-1/6); # sigma_OT in [A]
@p_water = (0,0,0, @p, -$p[0],$p[1],0);
$v_fcc = $n_water*$v_water; # cell volume
$l_fcc = $v_fcc**(1/3); # cell length
@p_fcc = (0.00,0.00,0.00, 0.50,0.50,0.00,
0.50,0.00,0.50, 0.00,0.50,0.50);
@n_fcc = ();
for (my $i=0; $i<scalar(@L); ++$i)
{
my $n = $L[$i]/$l_fcc; # calculate n_fcc
$n = int($n-int($n) ? $n+1 : $n); # ceil($n)
$L[$i] = $n*$l_fcc; # adjust box length
printf("Info: changed l%s to %g A\n", ("x","y","z")[$i], $L[$i])
if ($info);
push(@n_fcc, $n);
}
foreach (@p_fcc) { $_ = ($_+0.25)*$l_fcc; } # p_fcc in [A]
for (my $x=0; $x<$n_fcc[0]; ++$x) { # initialize flags
for (my $y=0; $y<$n_fcc[1]; ++$y) {
for (my $z=0; $z<$n_fcc[2]; ++$z) {
$flags_fcc[$x][$y][$z] = 15; } } } # turn on all fcc sites
}
sub floor
{
my $x = shift(@_);
return $x>0 ? int($x) : int($x)-1;
}
sub Periodic
{
my @p = splice(@_, 0, 3);
return (
$p[0]-floor($p[0]/$L[0]+0.5)*$L[0],
$p[1]-floor($p[1]/$L[1]+0.5)*$L[1],
$p[2]-floor($p[2]/$L[2]+0.5)*$L[2]);
}
sub EraseWater
{
my $r = shift(@_)/2;
my @p = splice(@_, 0, 3);
@p = V_Subtr(@p, @Center) if (!$center);
my @edges = (
$p[0]-$r,$p[1]-$r,$p[2]-$r, $p[0]-$r,$p[1]-$r,$p[2]+$r,
$p[0]-$r,$p[1]+$r,$p[2]-$r, $p[0]-$r,$p[1]+$r,$p[2]+$r,
$p[0]+$r,$p[1]-$r,$p[2]-$r, $p[0]+$r,$p[1]-$r,$p[2]+$r,
$p[0]+$r,$p[1]+$r,$p[2]-$r, $p[0]+$r,$p[1]+$r,$p[2]+$r);
my %list;
my @n;
my $d2 = ($r_water+$r)**2;
my @l = ($L[0]/2, $L[1]/2, $L[2]/2);
for (my $i=0; $i<scalar(@edges); $i+=3) # determine candidates
{
my @q = Periodic(@edges[$i, $i+1, $i+2]);
my @n = (int(($q[0]+$l[0])/$l_fcc),int(($q[1]+$l[1])/$l_fcc),
int(($q[2]+$l[2])/$l_fcc));
++$list{join(" ", @n)};
}
foreach (sort(keys(%list))) # check overlap
{
my @n = split(" ");
my @corner = ($n[0]*$l_fcc-$l[0]+$p_water[0],
$n[1]*$l_fcc-$l[1]+$p_water[1],
$n[2]*$l_fcc-$l[2]+$p_water[2]);
my $bit = 1;
my $flags = 0;
for (my $i=0; $i<scalar(@p_fcc); $i+=3)
{
my @q = V_Add(@corner, @p_fcc[$i,$i+1,$i+2]);
my @dp = Periodic(V_Subtr(@q, @p));
$flags |= $bit if (V_Dot(@dp, @dp)>$d2); # turn on fcc
$bit *= 2;
}
$flags_fcc[$n[0]][$n[1]][$n[2]] &= $flags; # set flags
}
}
sub CountFCC
{
my $n = 0;
return $n_fccs = 0 if (!$water_dens);
for (my $x=0; $x<$n_fcc[0]; ++$x) { # count water
for (my $y=0; $y<$n_fcc[1]; ++$y) {
for (my $z=0; $z<$n_fcc[2]; ++$z) {
my $bit = 1;
my $flags = $flags_fcc[$x][$y][$z];
for (my $i=0; $i<$n_water; ++$i) {
++$n if ($flags & $bit);
$bit *= 2; } } } }
return ($n_fccs = $n);
}
sub AddIons
{
my $n = ($n_waters = CountFCC())-int(abs($net_charge));
return if (!$ions);
printf("Warning: charge not neutralized: too little water\n") if ($n<0);
return if ($n<0);
printf(
"Warning: charge not neutralized: net charge (%g) is not an integer\n",
$net_charge) if ($net_charge!=int($net_charge));
my $n_na = $net_charge<0 ? int(abs($net_charge)) : 0;
my $n_cl = $net_charge>0 ? int(abs($net_charge)) : 0;
my $n_mol = int($ion_molar*$n*$v_water*1e-27*$nav+0.5);
my $n_atoms = ($n_na += $n_mol)+($n_cl += $n_mol);
$n_waters -= $n_atoms;
printf(
"Info: adding ions: [NaCl] = %g mol/l (%d Na+, %d Cl-)\n",
$n_mol/$n/$v_water/$nav/1e-27, $n_na, $n_cl) if ($info);
$n += int(abs($net_charge));
my $salt = 2**$n_water;
srand(time()); # seed random number
for (my $x=0; $x<$n_fcc[0]; ++$x) # replace water by ions
{
for (my $y=0; $y<$n_fcc[1]; ++$y)
{
for (my $z=0; $z<$n_fcc[2]; ++$z)
{
my $bit = 1;
my $flags = $flags_fcc[$x][$y][$z];
for (my $i=0; $i<$n_water; ++$i)
{
if ($flags & $bit)
{
my $prob = $n_atoms/$n;
--$n;
if (rand()<$prob)
{
my $na = rand()<$n_na/$n_atoms ? 1 : 0;
--$n_atoms;
if ($na) { --$n_na; } else { --$n_cl; }
$flags |= $salt*(1+$salt*$na)*$bit; # set type of ion
}
};
$bit *= 2;
}
$flags_fcc[$x][$y][$z] = $flags;
}
}
}
}
# LAMMPS output
sub WriteLAMMPSHeader # print lammps header
{
printf(LAMMPS "LAMMPS data file. %sCreated by $program v$version on %s\n",
($add ? "CGCMM Style. atom_style full. " : ""),`date`);
printf(LAMMPS "%12d atoms\n", $natoms);
printf(LAMMPS "%12d bonds\n", $nbonds);
printf(LAMMPS "%12d angles\n", $nangles);
printf(LAMMPS "%12d dihedrals\n", $ndihedrals);
printf(LAMMPS "%12d impropers\n\n", $nimpropers);
printf(LAMMPS "%12d atom types\n", $natom_types);
printf(LAMMPS "%12d bond types\n", $nbond_types);
printf(LAMMPS "%12d angle types\n", $nangle_types);
printf(LAMMPS "%12d dihedral types\n", $ndihedral_types);
printf(LAMMPS "%12d improper types\n\n", $nimproper_types);
}
sub WriteControlHeader
{
printf(PDB_CTRL "REMARK \n");
printf(PDB_CTRL "REMARK CONTROL PDB %s_ctrl.pdb\n", $project);
printf(PDB_CTRL "REMARK CREATED BY %s v%s ON %s",
$program, $version, `date`);
printf(PDB_CTRL "REMARK \n");
printf(PSF_CTRL "PSF\n");
printf(PSF_CTRL "\n");
printf(PSF_CTRL "%8d !NTITLE\n", 2);
printf(PSF_CTRL " REMARKS CONTROL PSF %s_ctrl.psf\n", $project);
printf(PSF_CTRL " REMARKS CREATED BY %s v%s ON %s",
$program, $version, `date`);
printf(PSF_CTRL "\n");
}
sub WriteBoxSize # print box limits
{
my @lo = V_Mult(@L[0,1,2], -1/2);
my @hi = V_Mult(@L[0,1,2], 1/2);
@lo = V_Add(@lo, @Center) if (!$center);
@hi = V_Add(@hi, @Center) if (!$center);
printf(LAMMPS "%12.8g %12.8g xlo xhi\n", $lo[0], $hi[0]);
printf(LAMMPS "%12.8g %12.8g ylo yhi\n", $lo[1], $hi[1]);
printf(LAMMPS "%12.8g %12.8g zlo zhi\n\n", $lo[2], $hi[2]);
}
sub WriteMasses # print mass list
{
my $k = 0;
printf(LAMMPS "Masses\n\n");
foreach (@types)
{
printf(LAMMPS "%8d %10.7g%s\n",
++$k, $masses{$_}, $add ? " # ".$ids{$_} : "");
}
printf(LAMMPS "\n");
}
sub WriteFCCAtoms
{
my $k = shift(@_);
my $res = shift(@_);
return $k if (!$water_dens);
$k_fcc = $k+1;
my @id = ($symbols{OT}, $symbols{HT}, $symbols{HT},
$symbols{SOD}, $symbols{CLA});
my @par = ();
my @charge = (-0.834, 0.417, 0.417, 1, -1);
my $salt = 2**$n_water;
my @l = ($L[0]/2, $L[1]/2, $L[2]/2);
my $iwater = 0;
my $isalt = 0;
foreach(@id) { push(@par, $link{$_}); }
for (my $x=0; $x<$n_fcc[0]; ++$x)
{
for (my $y=0; $y<$n_fcc[1]; ++$y)
{
for (my $z=0; $z<$n_fcc[2]; ++$z)
{
my @corner = ($x*$l_fcc-$l[0], $y*$l_fcc-$l[1], $z*$l_fcc-$l[2]);
my $flags = $flags_fcc[$x][$y][$z];
my $bit = 1;
for (my $i=0; $i<scalar(@p_fcc); $i+=3)
{
my $pair = $bit;
if ($flags & $pair)
{
my @p = V_Add(@corner, @p_fcc[$i,$i+1,$i+2]);
my $j = 0; # print water
my $n = scalar(@p_water);
++$res;
if ($flags & ($pair *= $salt)) # print salt ion
{ # sodium if highest
$j = $flags & ($pair*$salt) ? 3 : 4;
$n = 1;
$counter = ++$isalt;
}
else { $counter = ++$iwater; }
for (my $i=0; $i<$n; $i+=3)
{
my @xyz = V_Add(@p, @p_water[$i,$i+1,$i+2]);
@xyz = V_Add(@xyz, @Center) if (!$center);
printf(LAMMPS "%8d %7d %5d %9.6g %16.12g %16.12g %16.12g%s\n",
++$k, $res, $par[$j], $charge[$j], $xyz[0], $xyz[1],
$xyz[2], $add ? " # ".$types[$par[$j]-1] : "");
printf(PDB_CTRL "ATOM %6.6s %-4.4s %-3.3s %5.5s %3.3s ".
"%7.7s %7.7s %7.7s %5.5s %5.5s %4.4s %s\n", $k,
$types[$par[$j]-1], $n-1 ? "HOH" : "ION", $res, "",
$xyz[0], $xyz[1], $xyz[2], "1.00", "0.00", "",
$n-1 ? "WATR" : "SALT") if ($pdb_ctrl);
printf(PSF_CTRL "%8d %4.4s %-4.4s %-4.4s %-4.4s %4.4s ".
"%16.8e %7.7s %9.9s 0\n", $k, $n-1 ? "WATR" : "SALT",
$counter, $n-1 ? "HOH" : "ION", $types[$par[$j]-1], $id[$j],
$charge[$j], $masses{$id[$j]}, "") if ($pdb_ctrl);
++$j;
}
}
$bit *= 2;
}
}
}
}
return $k;
}
sub WritePSFAtoms()
{
my $n = PSFGoto(atoms);
my @res = (0, 0);
printf(PSF_CTRL "%8d !NATOM\n", $n+2*$n_waters+$n_fccs);
while (<PSF>)
{
last if (!$n--);
my @psf = split(" ");
if ($res[1]!=$psf[2]) { ++$res[0]; $res[1] = $psf[2]; }
printf(PSF_CTRL "%8d %4.4s %-4.4s %-4.4s %-4.4s %-4.4s ".
"%16.8e %7.7s %9.9s %s\n", $psf[0], $psf[1], $res[0],
$psf[3], $psf[4], $psf[5], $psf[6], $psf[7], "", $psf[8]);
}
}
sub WriteAtoms # print positions etc.
{
my $n = PSFGoto(atoms);
my $k = 0;
my @res = (0, 0);
CRDGoto(atoms);
$net_charge = 0;
printf(LAMMPS "Atoms%s\n\n",($add ? " # full" : "")) if ($n>0);
for (my $i=0; $i<$n; ++$i)
{
my @crd = $pdb ? NextPDB2CRD() : split(" ", <CRD>);
my @psf = split(" ", <PSF>);
my @xyz = MV_Dot(@R, @crd[-6, -5, -4]);
@xyz = V_Subtr(@xyz, @Center) if ($center);
if ($crd[-2]!=$res[1]) { ++$res[0]; $res[1] = $crd[-2]; }
printf(LAMMPS "%8d %7d %5d %9.6g %16.12g %16.12g %16.12g%s\n", ++$k,
$res[0], $link{$atom_types[$k]}, $psf[6], $xyz[0], $xyz[1], $xyz[2],
$add ? " # ".$types[$link{$atom_types[$k]}-1] : "");
printf(PDB_CTRL "ATOM %6.6s %-4.4s %-4.4s %4.4s %3.3s ".
"%7.7s %7.7s %7.7s %5.5s %5.5s %4.4s %s\n", $k,
$crd[-7], $crd[-8], $res[0], "", $xyz[0], $xyz[1], $xyz[2],
"1.00", $crd[-1], "", $crd[-3]) if ($pdb_ctrl);
next if (!$water_dens); # is water added?
$net_charge += $psf[6];
my @c = split(" ", $parms[$link{$atom_types[$k]}-1]);
EraseWater($c[1], @xyz);
}
$net_charge = int($net_charge*1e5+($net_charge>0?0.5:-0.5))/1e5;
AddIons() if ($water_dens);
WritePSFAtoms() if ($pdb_ctrl);
$k = WriteFCCAtoms($k, $res[0]+$res[1]);
printf(PDB_CTRL "END\n") if ($pdb_ctrl);
printf(LAMMPS "\n");
return $k;
}
sub WriteParameters # print parameters
{
my $mode = shift(@_)+1;
my $header = ("Pair","Bond","Angle","Dihedral","Improper")[$mode];
my $hint = ("# lj/charmm/coul/long", "# harmonic", "# charmm", "# charmm", "# harmonic")[$mode];
my $n = (4, 2, 4, 4, 2)[$mode];
my $k = 0;
printf("Info: converting ".lc($mode ? $header : "Atom")."s\n") if ($info);
if ($mode--)
{
BondedTypes($mode);
BondedParameters($mode);
%link = CrossLink(@types);
CorrectDihedralParameters() if ($mode==2);
@parms = Delete(1, @parms) if ($mode==3);
}
return 0 if (!scalar(@parms));
printf(LAMMPS "%s Coeffs %s\n\n", $header, ($add ? $hint : ""));
for (my $i=0; $i<scalar(@parms); ++$i)
{
if ($parms[$i] ne "")
{
foreach (split(":", $parms[$i]))
{
my @tmp = split(" ");
printf(LAMMPS "%8d", ++$k);
for (my $j=0; $j<$n; ++$j) {
printf(LAMMPS " %16.12g", $j<scalar(@tmp) ? $tmp[$j] : 0); }
printf(LAMMPS "%s\n", $add ? " # ".$types[$i] : "");
}
} else { ++$k; }
}
printf(LAMMPS "\n");
return $k;
}
sub WriteFCCBonded
{
my $mode = shift(@_);
my $k = shift(@_);
my $atom = $k_fcc;
return $k if (($mode>1)||!$water_dens);
my $type = $mode ? CreateID(HT, OT, HT) : CreateID(HT, OT);
my $id = $link{$type};
my $salt = 2**$n_water;
for (my $x=0; $x<$n_fcc[0]; ++$x)
{
for (my $y=0; $y<$n_fcc[1]; ++$y)
{
for (my $z=0; $z<$n_fcc[2]; ++$z)
{
my @corner = ($x*$l_fcc-$L[0]/2, $y*$l_fcc-$L[1]/2,
$z*$l_fcc-$L[2]/2);
my $flags = $flags_fcc[$x][$y][$z];
my $bit = 1;
for (my $i=0; $i<scalar(@p_fcc); $i+=3)
{
if ($flags&$bit)
{
if ($flags&($bit*$salt)) { ++$atom; }
else
{
printf(LAMMPS "%8d %7d %7d %7d%s\n", ++$k, $id, $atom,
$atom+1, $add ? " # ".$type : "") if (!$mode);
printf(LAMMPS "%8d %7d %7d %7d%s\n", ++$k, $id, $atom,
$atom+2, $add ? " # ".$type : "") if (!$mode);
printf(LAMMPS "%8d %7d %7d %7d %7d%s\n", ++$k, $id, $atom+1,
$atom, $atom+2, $add ? " # ".$type : "") if ($mode);
if ($pdb_ctrl)
{
PSFWrite(2, $atom, $atom+1, $atom, $atom+2) if (!$mode);
PSFWrite(3, $atom+1, $atom, $atom+2) if ($mode);
}
$atom += 3;
}
}
$bit *= 2;
}
}
}
}
return $k;
}
sub WriteBonded # print bonded list
{
my $mode = shift(@_);
my $psf_id = ("!NBOND:", "!NTHETA:", "!NPHI:", "!NIMPHI:")[$mode];
my $title = ("bonds", "angles", "dihedrals", "impropers")[$mode];
my $items = (2, 3, 4, 4)[$mode];
my $n = PSFGoto($title);
my $k = 0;
my @delta;
my @tmp;
return 0 if ($n<1);
printf(LAMMPS "%s\n\n", ucfirst($title));
printf(PSF_CTRL "\n%8d %s %s\n", $n+($mode ? ($mode==1 ? $n_waters : 0)
: 2*$n_waters), $psf_id, $title) if ($pdb_ctrl);
$psf_ncols = 0 if ($pdb_ctrl);
foreach (@parms)
{
push(@delta, $k);
$k += scalar(split(":"))-1 if ($_ ne "");
}
$k = 0;
for (my $i=0; $i<$n; ++$i)
{
my @bonded = PSFGet($items);
my @tmp = ();
foreach (@bonded) { push(@tmp, $ids{$atom_types[$_]}); }
my $id = $link{CreateID(@tmp)}-1;
my $m = 0;
if ($parms[$id] ne "")
{
foreach (split(":", $parms[$id]))
{
++$m;
my @const = split(" ");
next if (($const[0]==0)&&($mode==2 ? $const[-1]==0 : 1));
printf(LAMMPS "%8d %7d", ++$k, $id+$delta[$id]+$m);
foreach (@bonded) { printf(LAMMPS " %7d", $_); }
printf(LAMMPS "%s\n", $add ? " # ".CreateID(@tmp) : "");
}
}
else
{
printf(LAMMPS "%8d %7d", ++$k, $id+$delta[$id]+$m);
foreach (@bonded) { printf(LAMMPS " %7d", $_); }
printf(LAMMPS "%s\n", $add ? " # ".CreateID(@tmp) : "");
}
PSFWrite($items, @bonded) if ($pdb_ctrl);
}
$k = WriteFCCBonded($mode, $k);
printf(PSF_CTRL "\n") if ($pdb_ctrl && $psf_ncols);
printf(LAMMPS "\n");
return $k;
}
sub CreateCorrectedPairCoefficients
{
my $read = 0;
my $k = 0;
my %id;
my %type;
$coefficients = "";
foreach (@types) { $id{$ids{$_}} = $_; $type{$_} = ++$k; }
seek(PARAMETERS, 0, 0);
while (<PARAMETERS>)
{
chop();
my @cols = split(" ");
if ($read&&(scalar(@cols)>3)&&
(substr($cols[0],0,1) ne "!")&&($cols[2] lt 'A'))
{
my $id1 = $id{$cols[0]};
my $id2 = $id{$cols[1]};
if (($id1 ne "")&&($id2 ne ""))
{
my @c = (abs($cols[2]), $cols[3]*2.0**(-1/6));
if ($type{$id2}<$type{$id1})
{
my $tmp = $id1; $id1 = $id2; $id2 = $tmp;
}
$coefficients .= ":" if ($coefficients ne "");
$coefficients .= $type{$id1}." ".$type{$id2}." ";
$coefficients .= $c[0]." ".$c[1]." ".$c[0]." ".$c[1];
}
}
$read = 1 if ($cols[0] eq "NBFIX");
last if ($read&&!scalar(@cols));
}
}
sub WriteData
{
open(LAMMPS, ">$project.in"); # use .in for temporary
open(PDB_CTRL, ">".$project."_ctrl.pdb") if ($pdb_ctrl);
open(PSF_CTRL, ">".$project."_ctrl.psf") if ($pdb_ctrl);
WriteControlHeader() if ($pdb_ctrl);
ReadTopology();
CharacterizeBox();
SetupWater() if ($water_dens);
WriteBoxSize(); # body storage
@types = AtomTypes(); # atoms
@parms = AtomParameters(@types);
WriteMasses();
%link = CrossLink(@types);
CreateCorrectedPairCoefficients();
for (my $i=0; $i<scalar(@types); ++$i) { $types[$i] = $ids{$types[$i]}; }
$natom_types = WriteParameters(-1); # pairs
if ($#types+1 > $natom_types) {
print "Warning: $#types atom types present, but only $natom_types pair coeffs found\n";
# reset to what is found while determining the number of atom types.
$natom_types = $#types+1;
}
$natoms = WriteAtoms();
$nbond_types = WriteParameters(0); # bonds
$nbonds = WriteBonded(0);
$nangle_types = WriteParameters(1); # angles
$nangles = WriteBonded(1);
$shake = $link{CreateID(("HT", "OT", "HT"))};
$ndihedral_types = WriteParameters(2); # dihedrals
$ndihedrals = WriteBonded(2);
$nimproper_types = WriteParameters(3); # impropers
$nimpropers = WriteBonded(3);
close(LAMMPS); # close temp file
open(LAMMPS, ">$project.data"); # open data file
WriteLAMMPSHeader(); # header
open(TMP, "<$project.in"); # open temp file
while (<TMP>) { printf(LAMMPS "%s", $_); } # spool body
close(TMP); # close temp file
if ($pdb_ctrl)
{
#while (<PSF>) { printf(PSF_CTRL "%s", $_); }
close(PSF_CTRL); close(PDB_CTRL);
}
close(LAMMPS); # close data file
}
sub WriteLAMMPSInput
{
open(LAMMPS, ">$project.in"); # input file
printf(LAMMPS "# Created by $program v$version on %s", `date`);
printf(LAMMPS "# Command: %s\n\n", $cmdline);
printf(LAMMPS "units real\n"); # general
printf(LAMMPS "neigh_modify delay 2 every 1\n\n");
printf(LAMMPS "atom_style full\n"); # styles
printf(LAMMPS "bond_style harmonic\n") if ($nbond_types);
printf(LAMMPS "angle_style charmm\n") if ($nangle_types);
printf(LAMMPS "dihedral_style charmm\n") if ($ndihedral_types);
printf(LAMMPS "improper_style harmonic\n\n") if ($nimproper_types);
printf(LAMMPS "pair_style lj/charmm/coul/long 8 12\n");
printf(LAMMPS "pair_modify mix arithmetic\n");
printf(LAMMPS "kspace_style pppm 1e-6\n\n");
if ($cmap) {
printf(LAMMPS "# Modify following line to point to the desired CMAP file\n");
printf(LAMMPS "fix cmap all cmap charmm$cmap.cmap\n");
printf(LAMMPS "fix_modify cmap energy yes\n");
printf(LAMMPS "read_data $project.data fix cmap crossterm CMAP\n\n");
}else{
printf(LAMMPS "read_data $project.data\n\n"); # read data
}
if ($coefficients ne "") # corrected coeffs
{
foreach (split(":", $coefficients))
{
printf(LAMMPS "pair_coeff %s\n", $_);
}
printf(LAMMPS "\n");
}
printf(LAMMPS "special_bonds charmm\n"); # invoke charmm
printf(LAMMPS "thermo 10\n"); # set thermo style
printf(LAMMPS "thermo_style multi\n");
printf(LAMMPS "timestep 1.0\n\n"); # 1.0 ps time step
printf(LAMMPS "minimize 0.0 0.0 50 200\n\n"); # take of the edge
printf(LAMMPS "reset_timestep 0\n");
printf(LAMMPS "fix 1 all nve\n");
printf(LAMMPS "fix 2 all shake 1e-6 500 0 m 1.0\n")
if ($shake eq ""); # shake all H-bonds
printf(LAMMPS "fix 2 all shake 1e-6 500 0 m 1.0 a %s\n",$shake)
if ($shake ne ""); # add water if present
printf(LAMMPS "velocity all create 0.0 12345678 dist uniform\n\n");
printf(LAMMPS "restart 500 $project.restart1 $project.restart2\n");
printf(LAMMPS "dump 1 all atom 100 $project.dump\n");
printf(LAMMPS "dump_modify 1 image yes scale yes\n\n");
printf(LAMMPS "thermo 100\n"); # set thermo style
printf(LAMMPS "run 1000\n"); # run for 1000 time steps
close(LAMMPS);
}
# ----------------------- DESCRIPTION: sub CharmmCmap ------------------------ #
# This subroutine add a new section "CMAP" to the LAMMPS data file, which #
# a part of the implementation of "CHARMM CMAP" (see references) in LAMMPS. #
# The section "CMAP" contains a list of dihedral ID pairs from adjecent #
# peptide backtone dihedrals whose dihedral angles are corrresponding to PHI #
# and PSI. (PHI: C--N--C_aphla_C and PSI: N--C_alpha--C--N) #
# #
# Initiated by: Xiaohu Hu (hux2@ornl.gov) #
# May 2009 #
# #
# Finalized Oct 2016 by: Robert Latour (latourr@clemson.edu), #
# David Hyde-Volpe, and Tigran Abramyan, Clemson University, #
# and Chris Lorenz (chris.lorenz@kcl.ac.uk #
# #
# References: #
# - MacKerell, A.D., Jr., Feig, M., Brooks, C.L., III, Improved Treatment of #
# Protein Backbone Conformational in Empirical Force Fields, J. Am. Chem. #
# Soc. 126(2004): 698-699. #
# - MacKerell, A.D., Jr., Feig, M., Brooks, C.L., III, Extending the Treatment #
# of Backbone Energetics in Protein Force Fields: Limitations of Gas-Phase #
# Quantum Mechnacis in Reproducing Protein Conformational Distributions in #
# Molecular Dynamics Simulations, J. Comput. Chem. 25(2004): 1400-1415. #
# ---------------------------------------------------------------------------- #
sub CharmmCmap
{
print "\nINITIATING CHARMM CMAP SUBROUTINE...\n\n";
# Reread and analyse $project.data
my @raw_data;
open(LAMMPS, "< $project.data") or
die "\"sub CharmmCmap()\" cannot open \"$project.data!\n";
print "Analyzing \"$project.data\"...\n\n";
@raw_data = <LAMMPS>;
close(LAMMPS);
# Locate and extract the sections "Masses" and "Atoms"
my $line_number = 0;
# Header infos, 0 by default
my $natom_types = 0;
my $natom_number = 0;
my $ndihedral_number = 0;
my $temp_string;
# splice points, 0 by default
my $splice_onset_masses = 0;
my $splice_onset_atoms = 0;
my $splice_onset_dihedrals = 0;
foreach my $line (@raw_data) {
$line_number++;
chomp($line);
# Extract useful informations from the header
if ($line =~ m/atom types/) {
($natom_types,$temp_string) = split(" ",$line);
if ($natom_types == 0) {
die "\nError: Number of atom types is 0!\n";
}
print "Total atom types: $natom_types\n";
}
if ($line =~ m/atoms/) {
($natom_number,$temp_string) = split(" ",$line);
if ($natom_number == 0) {
die "\nError: Number of atoms is 0!\n";
}
print "Total atoms: $natom_number\n";
}
if ($line =~ m/dihedrals/) {
($ndihedral_number,$temp_string) = split(" ",$line);
if ($ndihedral_number == 0) {
die "\nError: Number of dihedrals is 0\n";
}
print "Total dihedrals: $ndihedral_number\n";
}
# Locate and data from sections "Masses", "Atoms" and "Dihedrals"
if ($line =~ m/Masses/) {
$splice_onset_masses = $line_number + 1;
if ($splice_onset_masses-1 == 0) {
die "\nError: Can not find the section \"Masses\"\n";
}
print "Section \"Masses\" found: line $splice_onset_masses\n";
}
if ($line =~ m/Atoms/) {
$splice_onset_atoms = $line_number +1;
if ($splice_onset_atoms-1 == 0) {
die "\nError: Can not find the section \"Atoms\"\n";
}
print "Section \"Atoms\" found: line $splice_onset_atoms\n";
}
if ($line =~ m/Dihedrals/) {
$splice_onset_dihedrals = $line_number + 1;
if ($splice_onset_dihedrals-1 == 0) {
die "\nError: Can not find the section \"Dihedrals\"\n";
}
print "Section \"Dihedrals\" found: line $splice_onset_dihedrals\n";
}
}
print "\nGenerating PHI/PSI dihedral pair list...\n\n";
my @temp1 = @raw_data;
my @temp2 = @raw_data;
my @temp3 = @raw_data;
# Extract the section "Masses", "Atoms" and "Dihedrals"
my @temp_masses_data = splice(@temp1,$splice_onset_masses,$natom_types);
my @temp_atoms_data = splice(@temp2,$splice_onset_atoms,$natom_number);
my @temp_dihedrals_data = splice(@temp3,$splice_onset_dihedrals,$ndihedral_number);
# Store @temp_masses_dat into a matrix
my @masses_matrix;
my $atom_type;
my $mass;
for (@temp_masses_data) {
($atom_type, $mass) = split(" ");
push(@masses_matrix,[$atom_type,$mass]);
}
# Store @temp_atoms_data into a matrix
my @atoms_matrix;
my $atom_ID;
my $molecule_ID;
my $atype;
my $charge;
my $atom_x_coor;
my $atom_y_coor;
my $atom_z_coor;
for (@temp_atoms_data) {
($atom_ID,$molecule_ID,$atype,$charge,$atom_x_coor,$atom_y_coor,$atom_z_coor) = split(" ");
push(@atoms_matrix,
[$atom_ID,$molecule_ID,$atype,$charge,$atom_x_coor,$atom_y_coor,$atom_z_coor]);
}
# Store @temp_dihedrals_data into a matrix
my @dihedrals_matrix;
my $dihedral_ID;
my $dihedtal_type;
my $dihe_atom1;
my $dihe_atom2;
my $dihe_atom3;
my $dihe_atom4;
for (@temp_dihedrals_data) {
($dihedral_ID,$dihedral_type,$dihe_atom1,$dihe_atom2,$dihe_atom3,$dihe_atom4) = split(" ");
push(@dihedrals_matrix,
[$dihedral_ID,$dihedral_type,$dihe_atom1,$dihe_atom2,$dihe_atom3,$dihe_atom4]);
}
# Find out and extract the peptide backbone dihedrals
#
# Definitions of peptide backbone dihedrals
#
# For dihedral angle PHI: C--N--CA--C
# For dihedral angle PSI: N--CA--C--N
#
# ---------------------------------------------------------
# atom | mass |partial charge| amino-acid
# ---------------------------------------------------------
# C | 12.011 | 0.51 | all except GLY and PRO
# N | 14.007 | -0.29 | PRO
# N | 14.007 | -0.47 | all except PRO
# CA | 12.011 | 0.07 | all except GLY and PRO
# CA | 12.011 | -0.02 | GLY
# CA | 12.011 | 0.02 | PRO
# ---------------------------------------------------------
#
# Peptide backbone
# ...
# /
# O=C
# \
# N-H
# / -----> PHI (C-N-CA-C)
# H-CA-R
# \ -----> PSI (N-CA-C-N)
# C=O
# /
# H-N
# \
# ...
#
# Criteria to be a PHI/PSI dihedral pair:
# 1. Atoms have to match with the mass/charge constellations as
# defined above.
# 2. The atoms N--CA--C needs to be covalently bonded with each
# other.
# Find which types do C, N and CA correspond to and store them
# in lists
my $mass_carbon = 12.011;
my $mass_nitrogen = 14.007;
my @carbon_list;
my @nitrogen_list;
my $carbon_counter = 0;
my $nitrogen_counter = 0;
for (my $i = 0; $i < $natom_types; $i++) {
if (${$masses_matrix[$i]}[1] == $mass_carbon) {
push(@carbon_list,${$masses_matrix[$i]}[0]);
$carbon_counter++;
}
if (${$masses_matrix[$i]}[1] == $mass_nitrogen) {
push(@nitrogen_list,${$masses_matrix[$i]}[0]);
$nitrogen_counter++;
}
}
# Quit if no carbons or nitrogens
if ($carbon_counter == 0 or $nitrogen_counter == 0) {
if ($carbon_counter == 0) {
print "No carbon atoms exist in the system\n";
}
if ($nitrogen_counter == 0) {
print "No nitrogen atoms exist in the system\n";
}
print "CMAP usage impossible\n";
return;
}
print "Carbon atom type/s: @carbon_list\n";
print "Nitrogen atom type/s: @nitrogen_list\n";
# Determine the atom types of C, CA and N
# Charges of the backbone atoms
my $charge_C = 0.51;
my $charge_CA = 0.07;
my $charge_N = -0.47;
# Special setting for PRO
my $charge_N_PRO = -0.29;
my $charge_CA_PRO = 0.02;
# Special setting for GLY
my $charge_CA_GLY = -0.02;
# Peptide backbone atom types
my $C_type;
my $CA_type;
my $CA_GLY_type;
my $CA_PRO_type;
my $N_type;
my $N_PRO_type;
my $C_counter = 0;
my $CA_counter = 0;
my $CA_GLY_counter = 0;
my $CA_PRO_counter = 0;
my $N_counter = 0;
my $N_PRO_counter = 0;
my $C_flag = 0;
for (my $i = 0; $i <= $natom_number; $i++) {
my $cur_type = ${$atoms_matrix[$i]}[2];
my $cur_charge = ${$atoms_matrix[$i]}[3];
for (my $j = 0; $j <= $#carbon_list; $j++) {
if ($cur_type == $carbon_list[$j]) {
$C_flag = 1;
if ($cur_charge == $charge_C) {
$C_type = $cur_type;
$C_counter++;
}
if ($cur_charge == $charge_CA) {
$CA_type = $cur_type;
$CA_counter++;
}
if ($cur_charge == $charge_CA_GLY) {
$CA_GLY_type = $cur_type;
$CA_GLY_counter++;
}
if ($cur_charge == $charge_CA_PRO) {
$CA_PRO_type = $cur_type;
$CA_PRO_counter++;
}
}
}
if ($C_flag == 0) {
for (my $k = 0; $k <= $#nitrogen_list; $k++) {
if ($cur_type == $nitrogen_list[$k]) {
if ($cur_charge == $charge_N) {
$N_type = $cur_type;
$N_counter++;
}
if ($cur_charge == $charge_N_PRO) {
$N_PRO_type = $cur_type;
$N_PRO_counter++;
}
}
}
}
$C_flag = 0;
}
# Quit if one of the atom types dosen't exist
if ( $C_counter == 0 or
($CA_counter == 0 and $CA_GLY_counter == 0 and $CA_PRO_counter == 0) or
($N_counter == 0 and $N_PRO_counter == 0) ) {
if ($C_counter == 0) {
print "\nCannot find the peptide backbone C atom type\n";
}
if ($CA_counter == 0 and $CA_GLY_counter == 0 and $CA_PRO_counter == 0) {
print "\nCannot find the peptide backbone C-alpha atom type\n";
}
if ($N_counter == 0 and $N_PRO_counter == 0) {
print "\nCannot find the peptide backbone N atom type\n";
}
print "CMAP usage impossible\n";
return;
}
print "Peptide backbone carbon type: $C_type\n";
print "Alpha-carbon type: $CA_type\n" if ($CA_counter > 0);
print "Alpha-carbon type (GLY): $CA_GLY_type\n" if ($CA_GLY_counter > 0);
print "Alpha-carbon type (PRO): $CA_PRO_type\n" if ($CA_PRO_counter > 0);
print "Peptide backbone nitrogen type: $N_type\n" if ($N_counter >0);
print "Peptide backbone nitrogen type (PRO): $N_PRO_type\n" if ($N_PRO_counter > 0);
# Loop through the dihedral list to find the PHI- and PSI-dihedrals
my @PHI_dihedrals;
my @PSI_dihedrals;
my $PHI_counter = 0;
my $PSI_counter = 0;
for (my $i = 0; $i < $ndihedral_number; $i++) {
my $cur_dihe_ID = ${dihedrals_matrix[$i]}[0];
my $cur_atom1_type = ${atoms_matrix[${dihedrals_matrix[$i]}[2]-1]}[2];
my $cur_atom2_type = ${atoms_matrix[${dihedrals_matrix[$i]}[3]-1]}[2];
my $cur_atom3_type = ${atoms_matrix[${dihedrals_matrix[$i]}[4]-1]}[2];
my $cur_atom4_type = ${atoms_matrix[${dihedrals_matrix[$i]}[5]-1]}[2];
next if (${dihedrals_matrix[$i]}[2] == ${dihedrals_matrix[$i-1]}[2] and
${dihedrals_matrix[$i]}[3] == ${dihedrals_matrix[$i-1]}[3] and
${dihedrals_matrix[$i]}[4] == ${dihedrals_matrix[$i-1]}[4] and
${dihedrals_matrix[$i]}[5] == ${dihedrals_matrix[$i-1]}[5]);
# Determine PHI-dihedrals; If C-CA-N-C or C-N-CA-C, then save it in a list
if ($cur_atom1_type == $C_type and $cur_atom4_type == $C_type) {
if ( ( ($cur_atom2_type == $CA_type or
$cur_atom2_type == $CA_GLY_type or
$cur_atom2_type == $CA_PRO_type) and
($cur_atom3_type == $N_type or
$cur_atom3_type == $N_PRO_type) ) or
( ($cur_atom3_type == $CA_type or
$cur_atom3_type == $CA_GLY_type or
$cur_atom3_type == $CA_PRO_type) and
($cur_atom2_type == $N_type or
$cur_atom2_type == $N_PRO_type) ) ) {
push (@PHI_dihedrals,$cur_dihe_ID);
$PHI_counter++;
}
}
# Determin PSI-dihedrals; If N-CA-C-N or N-C-CA-N (N can be both normal N or N proline),
# then save it in a list
if ( ($cur_atom1_type == $N_type and $cur_atom4_type == $N_type) or
($cur_atom4_type == $N_PRO_type and $cur_atom1_type == $N_PRO_type) or
($cur_atom1_type == $N_type and $cur_atom4_type == $N_PRO_type) or
($cur_atom4_type == $N_type and $cur_atom1_type == $N_PRO_type) ) {
if ( ( ($cur_atom2_type == $CA_type or
$cur_atom2_type == $CA_GLY_type or
$cur_atom2_type == $CA_PRO_type) and
$cur_atom3_type == $C_type) or
( ($cur_atom3_type == $CA_type or
$cur_atom3_type == $CA_GLY_type or
$cur_atom3_type == $CA_PRO_type) and
$cur_atom2_type == $C_type) ) {
push (@PSI_dihedrals,$cur_dihe_ID);
$PSI_counter++;
}
}
}
# Quit if no PHI or PSI dihedrals
if ($PHI_counter == 0 or $PSI_counter ==0) {
if ($PHI_counter == 0) {
print "Can not find the PHI backbone dihedrals\n";
}
if ($PSI_counter ==0) {
print "Can not find the PSI backbone dihedrals\n";
}
print "CMAP usage impossible\n";
return;
}
# Construct the PHI/PSI diheral pair list
#
# The algorithm:
# _____
# | |
# 1--2--3--4 PHI-dihedral
# 4--3--2--1
# --C--N-CA--C--N-- Peptide backbone
# 1--2--3--4
# 4--3--2--1 PSI-dihedral
# |_____|
#
# For a certain PHI dihedral, following conditions have to be met:
#
# PHI PSI
# If (2--3--4) = (1--2--3)
# or
# if (2--3--4) = (4--3--2)
# or
# if (3--2--1) = (1--2--3)
# or
# if (3--2--1) = (4--3--2),
#
# then these 2 dihedrals are a PHI/PSI pair. If a pair is found, the
# dihedral IDs will be stored in "@PHI_PSI_matrix".
my @PHI_PSI_matrix;
my $crossterm_CA_charge;
my $crossterm_type;
my $crossterm_counter = 0;
my $crossterm_type1_flag = 0;
my $crossterm_type2_flag = 0;
my $crossterm_type3_flag = 0;
my $crossterm_type4_flag = 0;
my $crossterm_type5_flag = 0;
my $crossterm_type6_flag = 0;
for (my $i = 0; $i <= $#PHI_dihedrals; $i++) {
my $cur_PHI_dihe = ${dihedrals_matrix[$PHI_dihedrals[$i]-1]}[0];
my $phi1 = ${dihedrals_matrix[$PHI_dihedrals[$i]-1]}[2];
my $phi2 = ${dihedrals_matrix[$PHI_dihedrals[$i]-1]}[3];
my $phi3 = ${dihedrals_matrix[$PHI_dihedrals[$i]-1]}[4];
my $phi4 = ${dihedrals_matrix[$PHI_dihedrals[$i]-1]}[5];
for (my $j = 0; $j <= $#PSI_dihedrals; $j++) {
my $cur_PSI_dihe = ${dihedrals_matrix[$PSI_dihedrals[$j]-1]}[0];
my $psi1 = ${dihedrals_matrix[$PSI_dihedrals[$j]-1]}[2];
my $psi2 = ${dihedrals_matrix[$PSI_dihedrals[$j]-1]}[3];
my $psi3 = ${dihedrals_matrix[$PSI_dihedrals[$j]-1]}[4];
my $psi4 = ${dihedrals_matrix[$PSI_dihedrals[$j]-1]}[5];
if ( ($phi2 == $psi1 and $phi3 == $psi2 and $phi4 == $psi3) or
($phi2 == $psi4 and $phi3 == $psi3 and $phi4 == $psi2) or
($phi3 == $psi1 and $phi2 == $psi2 and $phi1 == $psi3) or
($phi3 == $psi4 and $phi2 == $psi3 and $phi1 == $psi2) ) {
# Find out to which amino acid the cross-term belongs
if ($phi3 == $psi2 or $phi3 == $psi3) {
$crossterm_CA_charge = ${atoms_matrix[${dihedrals_matrix[$PHI_dihedrals[$i]-1]}[4]-1]}[3];
}
if ($phi2 == $psi2 or $phi2 == $psi3) {
$crossterm_CA_charge = ${atoms_matrix[${dihedrals_matrix[$PHI_dihedrals[$i]-1]}[3]-1]}[3];
}
# Def. the type of the crossterm per cmap.data file; If C_alpha of the crossterm is
# - ALA type, then $crossterm_type = 1;
# - ALA-PRO (ALA is the current AA), then $crossterm_type = 2;
# - PRO type, then $crossterm_type = 3;
# - PRO-PRO (First PRO is the current AA), then $crossterm_type = 4;
# - GLY type, then $crossterm_type = 5;
# - GLY-PRO (GLY is the current AA), then $crossterm_type = 6;
if ($crossterm_CA_charge == $charge_CA) { $crossterm_type = 1; $crossterm_type1_flag = 1; }
if ($crossterm_CA_charge == $charge_CA_GLY) { $crossterm_type = 5; $crossterm_type5_flag = 1; }
if ($crossterm_CA_charge == $charge_CA_PRO) {
$crossterm_type = 3; $crossterm_type3_flag = 1;
# Checking the last crossterm, re-assign the last crossterm type if needed
if ($crossterm_counter-1 >= 0 and $PHI_PSI_matrix[$crossterm_counter-1][0] == 1) {
$PHI_PSI_matrix[$crossterm_counter-1][0] = 2;
$crossterm_type2_flag = 1;
}
if ($crossterm_counter-1 >= 0 and $PHI_PSI_matrix[$crossterm_counter-1][0] == 3) {
$PHI_PSI_matrix[$crossterm_counter-1][0] = 4;
$crossterm_type4_flag = 1;
}
if ($crossterm_counter-1 >= 0 and $PHI_PSI_matrix[$crossterm_counter-1][0] == 5) {
$PHI_PSI_matrix[$crossterm_counter-1][0] = 6;
$crossterm_type6_flag = 1;
}
}
push(@PHI_PSI_matrix,[$crossterm_type,$phi1,$phi2,$phi3,$phi4,$psi4]);
$crossterm_counter++;
$crossterm_CA_charge = 0;
$crossterm_type = 0;
}
}
}
# Check whether the amino acid at the C-terminus is a PRO or not. If yes, the type of the last crossterm
# should be set to its X-PRO form instead of X, where X is ALA, PRO, or GLY. X-PRO form = X type + 1.
my @pdb_data;
open(PDB,"< $project.pdb")
or die "WARNING: Cannot open file \"$project.pdb\"! (required if the -cmap option is used)\n";
@pdb_data = <PDB>;
close(PDB);
my @ter_line;
my $ter_AA_type = 0;
my $ter_flag = 0;
foreach $line (@pdb_data) {
if ($line =~ m/TER/) {
@ter_line = split(" ",$line);
$ter_AA_type = $ter_line[2];
print "Terminal amino acid type is: $ter_AA_type\n";
$ter_flag = 1;
}
}
if ($ter_flag == 0) {
print "\n*** ERROR IN THE PDB FILE: ***\n";
print "In order for the CMAP section to be generated, the pdb file must \n";
print "identify the C-terminus amino acid in the file with 'TER'. \n";
print "This line is missing from the pdb file that was used.\n";
print "To correct this problem, open the pdb file in an editor,\n";
print "find the last atom of the last amino acid residue in the peptide\n";
print "chain and insert the following line immediately after that atom:\n";
print " 'TER <#1> <RES> <#2>' \n";
print "where '<#1> is the next atom number, <RES> is the three letter amino\n";
print "acid abbreviation for that amino acid, and <#2> is the molecule number\n";
print "of the terminal amino acid residue.\n\n";
print "For example, if the last atom of the last amino acid in the peptide\n";
print "sequence is listed in the pdb file as:\n\n";
print " 'ATOM 853 O GLU P 56 12.089 -1.695 -6.543 1.00 1.03 PROA'\n\n";
print "you would insert the following line after it:\n\n";
print " 'TER 854 GLU 56'\n\n";
print "If any additional atoms are listed in the pdb file (e.g., water, ions)\n";
print "after this terminal amino acid residue, their atom numbers and\n";
print "molecule numbers must be incremented by 1 to account for the new line\n";
print "that was inserted.\n\n";
die "Error: No terminating atom designated in pdb file! See above note to correct problem.\n\n";
}
if ($ter_AA_type eq PRO) {
$PHI_PSI_matrix[$crossterm_counter-1][0] = $PHI_PSI_matrix[$crossterm_counter-1][0]+1;
}
# Print out PHI/PSI diheral pair list
my $pair_counter = 0;
# Don't presently use $ncrosstermtypes but have this available if wish to print it out
my $ncrosstermtypes = $crossterm_type1_flag + $crossterm_type2_flag + $crossterm_type3_flag +
$crossterm_type4_flag + $crossterm_type5_flag + $crossterm_type6_flag;
print "\nWriting \"$project.data\" with section \"CMAP crossterms\" added at the end.\n";
# Writing the new lammps data file
open(REWRITE,"> $project.data")
or die "Cannot write file \"$project.data\"!\n";
foreach $line (@raw_data) {
printf(REWRITE "$line\n");
if ($line =~ m/impropers/) {
printf(REWRITE "%12d crossterms\n", $crossterm_counter);
}
}
printf(REWRITE "CMAP\n\n");
my $ref_line;
my $column;
foreach $ref_line (@PHI_PSI_matrix) {
$pair_counter++;
printf(REWRITE "%8d",$pair_counter);
foreach $column (@$ref_line) {
printf(REWRITE " %7d",$column);
}
printf(REWRITE "\n");
}
close(REWRITE);
print "\nDone!\n\n";
}
# main
Initialize();
WriteData();
WriteLAMMPSInput();
printf("Info: conversion complete\n\n") if ($info);
CharmmCmap() if ($cmap);
Event Timeline
Log In to Comment