diff --git a/tools/ch2lmp/charmm2lammps.pl b/tools/ch2lmp/charmm2lammps.pl
index dbc6ba003cd042e485122aa1455b2e863e9c01aa..17b4b817aff5d8f058727e39a33eff779817a8f7 100755
--- a/tools/ch2lmp/charmm2lammps.pl
+++ b/tools/ch2lmp/charmm2lammps.pl
@@ -1,33 +1,33 @@
#!/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
+# 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
+# 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
@@ -35,70 +35,71 @@
# 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
+# 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".
+# 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(@_);
+ my $name = shift(@_);
printf("Error: file %s not found\n", $name) if (!scalar(stat($name)));
return !scalar(stat($name));
}
- sub Initialize # initialization
+ 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",
+ 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"
- );
+ "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.0";
- $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();
+ $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";
@@ -120,45 +121,48 @@
$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
+ 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);
+ $forcefield = $_ if (!$k);
+ $project = $_ if ($k++ == 1);
}
}
- $water_dens = 1 if ($ions && !$water_dens);
+ $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",
@@ -167,25 +171,25 @@
printf("Options:\n");
for (my $i=0; $i<scalar(@options); ++$i)
{
- printf(" %-10.10s %s\n", $options[$i], $remarks[$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;
+ 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]);
+ ("x", "y", "z")[$i]) if ($info&&!$L[$i]);
}
open(PARAMETERS, "<par_$forcefield.prm");
}
@@ -195,7 +199,7 @@
sub V_String
{
- my @v = @_;
+ my @v = @_;
return "{".$v[0].", ".$v[1].", ".$v[2]."}";
}
@@ -203,8 +207,8 @@
sub V_Add
{
- my @v1 = splice(@_, 0, 3);
- my @v2 = splice(@_, 0, 3);
+ my @v1 = splice(@_, 0, 3);
+ my @v2 = splice(@_, 0, 3);
return ($v1[0]+$v2[0], $v1[1]+$v2[1], $v1[2]+$v2[2]);
}
@@ -212,8 +216,8 @@
sub V_Subtr
{
- my @v1 = splice(@_, 0, 3);
- my @v2 = splice(@_, 0, 3);
+ my @v1 = splice(@_, 0, 3);
+ my @v2 = splice(@_, 0, 3);
return ($v1[0]-$v2[0], $v1[1]-$v2[1], $v1[2]-$v2[2]);
}
@@ -221,8 +225,8 @@
sub V_Dot
{
- my @v1 = splice(@_, 0, 3);
- my @v2 = splice(@_, 0, 3);
+ my @v1 = splice(@_, 0, 3);
+ my @v2 = splice(@_, 0, 3);
return $v1[0]*$v2[0]+$v1[1]*$v2[1]+$v1[2]*$v2[2];
}
@@ -230,8 +234,8 @@
sub V_Mult
{
- my @v = splice(@_, 0, 3);
- my $f = shift(@_);
+ my @v = splice(@_, 0, 3);
+ my $f = shift(@_);
return ($f*$v[0], $f*$v[1], $f*$v[2]);
}
@@ -243,8 +247,8 @@
for (my $i=0; $i<3; ++$i)
{
- $string .= ", " if ($i);
- $string .= V_String(splice(@_, 0, 3));
+ $string .= ", " if ($i);
+ $string .= V_String(splice(@_, 0, 3));
}
return "{".$string."}";
}
@@ -261,13 +265,13 @@
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);
+ 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),
@@ -281,14 +285,14 @@
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);
+ { # 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
@@ -298,10 +302,10 @@
sub MV_Dot
{
- my @v11 = splice(@_, 0, 3);
- my @v12 = splice(@_, 0, 3);
- my @v13 = splice(@_, 0, 3);
- my @v2 = splice(@_, 0, 3);
+ 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));
}
@@ -311,55 +315,55 @@
sub PSFConnectivity
{
- my $n = PSFGoto(bonds);
+ 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);
+ 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
+ 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 = ();
+ my $n = scalar(@nconnect);
+ my @bonded = ();
for (my $i=1; $i<=$n; ++$i)
{
- $bonded[0] = $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);
- }
- }
+ $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;
+ $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;
+ 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);
@@ -368,20 +372,20 @@
while (<PSF>)
{
chop();
- my @tmp = split(" ");
- my $n = $hash{$tmp[1]};
- $PSFIndex{$n} = tell(PSF)." ".$tmp[0] if ($n ne "");
+ my @tmp = split(" ");
+ my $n = $hash{$tmp[1]};
+ $PSFIndex{$n} = tell(PSF)." ".$tmp[0] if ($n ne "");
}
}
- sub PSFGoto # goto $ident in <PSF>
+ sub PSFGoto # goto $ident in <PSF>
{
CreatePSFIndex() if (!scalar(%PSFIndex));
- my $id = shift(@_);
- my @n = split(" ", $PSFIndex{$id});
+ my $id = shift(@_);
+ my @n = split(" ", $PSFIndex{$id});
- @PSFBuffer = ();
+ @PSFBuffer = ();
# return PSFDihedrals() if ($id eq "dihedrals");
if (!scalar(@n))
{
@@ -398,14 +402,14 @@
{
if ($dihedral_flag)
{
- $dihedral_flag = $idihedral+1<$ndihedral ? 1 : 0;
+ $dihedral_flag = $idihedral+1<$ndihedral ? 1 : 0;
return split(" ", $dihedral[$idihedral++]);
}
if (!scalar(@PSFBuffer))
{
- my $line = <PSF>;
+ my $line = <PSF>;
chop($line);
- @PSFBuffer = split(" ", $line);
+ @PSFBuffer = split(" ", $line);
}
return splice(@PSFBuffer, 0, shift(@_));
}
@@ -413,8 +417,8 @@
sub PSFWrite
{
- my $items = shift(@_);
- my $n = $items;
+ my $items = shift(@_);
+ my $n = $items;
if ($psf_ncols>7) { printf(PSF_CTRL "\n"); $psf_ncols = 0; }
foreach(@_)
@@ -423,9 +427,9 @@
++$psf_ncols;
if ((!--$n) && ($psf_ncols>7))
{
- printf(PSF_CTRL "\n");
- $psf_ncols = 0;
- $n = $items;
+ printf(PSF_CTRL "\n");
+ $psf_ncols = 0;
+ $n = $items;
}
}
}
@@ -461,36 +465,36 @@
sub Delete
{
- my $item = shift(@_);
- my $k = 0;
+ my $item = shift(@_);
+ my $k = 0;
my @list;
foreach (@_)
{
- my @tmp = split(" ");
+ my @tmp = split(" ");
delete($tmp[$item]);
- $list[$k++] = join(" ", @tmp);
+ $list[$k++] = join(" ", @tmp);
}
return @list;
}
- sub CreateID # create id from 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 $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]));
+ $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));
+ $id .= ($i ? " " : "").($tmp = $list[$flag ? --$j : $i]);
+ $id .= substr(" ", 0, 4-length($tmp));
}
return $id;
}
@@ -498,15 +502,15 @@
sub AtomTypes
{
- my $n = PSFGoto(atoms);
+ my $n = PSFGoto(atoms);
my %list;
return () if ($n<1);
- $atom_types[0] = -1;
+ $atom_types[0] = -1;
for (my $i=0; $i<$n; ++$i)
{
- my @tmp = split(" ", <PSF>);
- $tmp[5] = $symbols{$tmp[5]}
+ 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]};
@@ -541,44 +545,44 @@
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
+ 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
+ sub AtomParameters # non-bonded parameters
{
my @types;
my @list;
- my $k = 0;
- my $read = 0;
- my %markers = Markers();
+ my $k = 0;
+ my $read = 0;
+ my %markers = Markers();
foreach(@_) { $types{$ids{$_}} = $k++; }
seek(PARAMETERS, 0, 0);
while (<PARAMETERS>)
{
chop();
- my @cols = split(" ");
+ my @cols = split(" ");
if ($read&&(scalar(@cols)>1)&&
- (substr($cols[0],0,1) ne "!")&&($cols[1] lt "A"))
+ (substr($cols[0],0,1) ne "!")&&($cols[1] lt "A"))
{
- my $k = $types{shift(@cols)};
- $list[$k] = NonBond(@cols) if ($k ne "");
+ 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; }
+ $read = ($markers{$cols[0]} eq "0") ? 1 : 0; }
}
- $list[$types{HT}] = NonBond(0, -0.046, 0.2245)
+ $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)
+ $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)
+ $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 ""));
@@ -586,103 +590,103 @@
}
- sub BondedTypes # create bonded types
+ 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 $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 = ();
+ 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));
+ @types = sort(keys(%list));
}
- sub Parameters # parms from columns
+ sub Parameters # parms from columns
{
- my $items = shift(@_);
- my @cols = @_;
- my $parms = "";
+ 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];
+ $parms = $parms.($i>$items ? " " : "").$cols[$i];
}
return $parms;
}
- sub BondedParameters # distil parms from
- { # <PARAMETERS>
- my $mode = shift(@_); # bonded mode
+ 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 $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 = ();
+ @parms = ();
foreach(@types) { $link{$_} = $k++; }
seek(PARAMETERS, 0, 0);
while (<PARAMETERS>)
{
chomp();
- my @cols = split(" ");
+ 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 (($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; }
+ $read = ($markers{$cols[0]} eq $mode+1) ? 1 : 0; }
}
if ($water_dens)
{
@@ -692,31 +696,31 @@
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 "");
+ $name, $i+1, $types[$i]) if ($parms[$i] eq "");
}
}
- sub SetScreeningFactor # set screening factor
+ sub SetScreeningFactor # set screening factor
{
- my $id = shift(@_);
- my $value = shift(@_);
- my $new = "";
+ 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);
+ my @tmp = split(" ");
+ $tmp[-1] = $value if ($tmp[-1]);
+ $new .= ":" if ($new ne "");
+ $new .= join(" ", @tmp);
}
- $parms[$id] = $new;
+ $parms[$id] = $new;
}
sub CorrectDihedralParameters
{
- my $n = PSFGoto(dihedrals);
+ my $n = PSFGoto(dihedrals);
my %hash;
my $hash_id;
my $id1;
@@ -726,33 +730,33 @@
for (my $i=0; $i<$n; ++$i)
{
- my @bonded = PSFGet(4);
- my @tmp = ();
+ my @bonded = PSFGet(4);
+ my @tmp = ();
foreach (@bonded) { push(@tmp, $ids{$atom_types[$_]}); }
- $id1 = $link{CreateID(@tmp)}-1;
- $first = $bonded[0];
- $last = $bonded[3];
+ $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
+ $hash{$hash_id} = $id1; # add id to hash
}
else
{
- SetScreeningFactor($id1, 0.5); # 6-ring: shared 1-4
- SetScreeningFactor($id2, 0.5);
+ SetScreeningFactor($id1, 0.5); # 6-ring: shared 1-4
+ SetScreeningFactor($id2, 0.5);
}
}
- $n = PSFGoto(angles);
+ $n = PSFGoto(angles);
for (my $i=0; $i<$n; ++$i)
{
- my @bonded = PSFGet(3);
- $first = $bonded[0];
- $last = $bonded[2];
+ 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
+ SetScreeningFactor($id1, 0); # 5-ring: no 1-4
}
}
}
@@ -760,39 +764,39 @@
sub AddMass
{
- my $symbol = shift(@_);
- my $mass = shift(@_);
+ my $symbol = shift(@_);
+ my $mass = shift(@_);
return if ($symbols{$symbol} ne "");
- $ids{++$max_id} = $symbol;
- $masses{$max_id} = $mass;
- $symbols{$symbol} = $max_id;
+ $ids{++$max_id} = $symbol;
+ $masses{$max_id} = $mass;
+ $symbols{$symbol} = $max_id;
}
- sub ReadTopology # read topology links
+ sub ReadTopology # read topology links
{
- my $id = shift(@_);
- my $item = shift(@_);
- my $read = 0;
+ my $id = shift(@_);
+ my $item = shift(@_);
+ my $read = 0;
my @tmp;
open(TOPOLOGY, "<top_$forcefield.rtf");
- $max_id = 0;
+ $max_id = 0;
while (<TOPOLOGY>)
{
- chop(); # delete CR at end
- my @tmp = split(" ");
- $read = 1 if ($tmp[0] eq "MASS");
+ 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]);
+ $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
+ # $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);
@@ -802,10 +806,10 @@
}
- sub CrossLink # symbolic cross-links
+ sub CrossLink # symbolic cross-links
{
- my @list = @_;
- my $n = scalar(@list);
+ my @list = @_;
+ my $n = scalar(@list);
my %hash;
for (my $i=0; $i<$n; ++$i) { $hash{$list[$i]} = $i+1; }
@@ -815,35 +819,35 @@
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);
+ 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);
+ 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);
+ $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);
}
@@ -852,44 +856,44 @@
{
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);
+ 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 = ();
+ $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
+ 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);
+ 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
+ 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
+ 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(@_);
+ my $x = shift(@_);
return $x>0 ? int($x) : int($x)-1;
}
@@ -897,7 +901,7 @@
sub Periodic
{
- my @p = splice(@_, 0, 3);
+ my @p = splice(@_, 0, 3);
return (
$p[0]-floor($p[0]/$L[0]+0.5)*$L[0],
@@ -908,10 +912,10 @@
sub EraseWater
{
- my $r = shift(@_)/2;
- my @p = splice(@_, 0, 3);
- @p = V_Subtr(@p, @Center) if (!$center);
- my @edges = (
+ 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,
@@ -919,55 +923,55 @@
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 $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));
+ 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
+ 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;
+ 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;
+ 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
+ $flags_fcc[$n[0]][$n[1]][$n[2]] &= $flags; # set flags
}
}
sub CountFCC
{
- my $n = 0;
+ my $n = 0;
return $n_fccs = 0 if (!$water_dens);
- for (my $x=0; $x<$n_fcc[0]; ++$x) { # count water
+ 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; } } } }
+ 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));
+ my $n = ($n_waters = CountFCC())-int(abs($net_charge));
return if (!$ions);
printf("Warning: charge not neutralized: too little water\n") if ($n<0);
@@ -975,43 +979,43 @@
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;
+ 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
+ $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;
- }
+ 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;
+ }
}
}
}
@@ -1019,7 +1023,7 @@
# LAMMPS output
- sub WriteLAMMPSHeader # print lammps header
+ 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`);
@@ -1054,28 +1058,28 @@
}
- sub WriteBoxSize # print box limits
+ 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);
+ 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);
+ @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
+ sub WriteMasses # print mass list
{
- my $k = 0;
+ my $k = 0;
printf(LAMMPS "Masses\n\n");
foreach (@types)
{
printf(LAMMPS "%8d %10.7g%s\n",
- ++$k, $masses{$_}, $add ? " # ".$ids{$_} : "");
+ ++$k, $masses{$_}, $add ? " # ".$ids{$_} : "");
}
printf(LAMMPS "\n");
}
@@ -1083,67 +1087,67 @@
sub WriteFCCAtoms
{
- my $k = shift(@_);
- my $res = shift(@_);
+ 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;
+ $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;
- }
- }
+ 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;
@@ -1152,76 +1156,76 @@
sub WritePSFAtoms()
{
- my $n = PSFGoto(atoms);
- my @res = (0, 0);
+ 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(" ");
+ 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]);
+ "%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.
+ sub WriteAtoms # print positions etc.
{
- my $n = PSFGoto(atoms);
- my $k = 0;
- my @res = (0, 0);
+ my $n = PSFGoto(atoms);
+ my $k = 0;
+ my @res = (0, 0);
CRDGoto(atoms);
- $net_charge = 0;
+ $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);
+ 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] : "");
+ $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]);
+ "%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;
+ $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]);
+ $k = WriteFCCAtoms($k, $res[0]+$res[1]);
printf(PDB_CTRL "END\n") if ($pdb_ctrl);
printf(LAMMPS "\n");
return $k;
}
- sub WriteParameters # print parameters
+ sub WriteParameters # print parameters
{
- my $mode = shift(@_)+1;
- my $header = ("Pair","Bond","Angle","Dihedral","Improper")[$mode];
+ 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;
+ 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);
+ %link = CrossLink(@types);
CorrectDihedralParameters() if ($mode==2);
- @parms = Delete(1, @parms) if ($mode==3);
+ @parms = Delete(1, @parms) if ($mode==3);
}
return 0 if (!scalar(@parms));
printf(LAMMPS "%s Coeffs %s\n\n", $header, ($add ? $hint : ""));
@@ -1229,14 +1233,14 @@
{
if ($parms[$i] ne "")
{
- foreach (split(":", $parms[$i]))
- {
- my @tmp = split(" ");
+ 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 " %16.12g", $j<scalar(@tmp) ? $tmp[$j] : 0); }
printf(LAMMPS "%s\n", $add ? " # ".$types[$i] : "");
- }
+ }
} else { ++$k; }
}
printf(LAMMPS "\n");
@@ -1246,104 +1250,104 @@
sub WriteFCCBonded
{
- my $mode = shift(@_);
- my $k = shift(@_);
- my $atom = $k_fcc;
+ 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;
+ 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;
- }
- }
+ 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
+ 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 $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);
+ : 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 += scalar(split(":"))-1 if ($_ ne "");
}
- $k = 0;
+ $k = 0;
for (my $i=0; $i<$n; ++$i)
{
- my @bonded = PSFGet($items);
- my @tmp = ();
+ my @bonded = PSFGet($items);
+ my @tmp = ();
foreach (@bonded) { push(@tmp, $ids{$atom_types[$_]}); }
- my $id = $link{CreateID(@tmp)}-1;
- my $m = 0;
+ 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));
+ ++$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", $_); }
+ 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", $_); }
+ foreach (@bonded) { printf(LAMMPS " %7d", $_); }
printf(LAMMPS "%s\n", $add ? " # ".CreateID(@tmp) : "");
}
PSFWrite($items, @bonded) if ($pdb_ctrl);
}
- $k = WriteFCCBonded($mode, $k);
+ $k = WriteFCCBonded($mode, $k);
printf(PSF_CTRL "\n") if ($pdb_ctrl && $psf_ncols);
printf(LAMMPS "\n");
return $k;
@@ -1352,36 +1356,36 @@
sub CreateCorrectedPairCoefficients
{
- my $read = 0;
- my $k = 0;
+ my $read = 0;
+ my $k = 0;
my %id;
my %type;
- $coefficients = "";
+ $coefficients = "";
foreach (@types) { $id{$ids{$_}} = $_; $type{$_} = ++$k; }
seek(PARAMETERS, 0, 0);
while (<PARAMETERS>)
{
chop();
- my @cols = split(" ");
+ my @cols = split(" ");
if ($read&&(scalar(@cols)>3)&&
- (substr($cols[0],0,1) ne "!")&&($cols[2] lt 'A'))
+ (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];
- }
+ 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");
+ $read = 1 if ($cols[0] eq "NBFIX");
last if ($read&&!scalar(@cols));
}
}
@@ -1389,58 +1393,59 @@
sub WriteData
{
- open(LAMMPS, ">$project.in"); # use .in for temporary
+ 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);
+ WriteBoxSize(); # body storage
+ @types = AtomTypes(); # atoms
+ @parms = AtomParameters(@types);
WriteMasses();
- %link = CrossLink(@types);
+ %link = CrossLink(@types);
CreateCorrectedPairCoefficients();
for (my $i=0; $i<scalar(@types); ++$i) { $types[$i] = $ids{$types[$i]}; }
- $natom_types = WriteParameters(-1); # pairs
+ $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
+ $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
+ close(LAMMPS); # close data file
}
sub WriteLAMMPSInput
{
- open(LAMMPS, ">$project.in"); # input file
- printf(LAMMPS "# Created by $program v$version on %s\n", `date`);
- printf(LAMMPS "units real\n"); # general
+ 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 "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);
@@ -1455,32 +1460,34 @@
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
+ printf(LAMMPS "read_data $project.data\n\n"); # read data
}
- if ($coefficients ne "") # corrected coeffs
+ if ($coefficients ne "") # corrected coeffs
{
foreach (split(":", $coefficients))
{
- printf(LAMMPS "pair_coeff %s\n", $_);
+ printf(LAMMPS "pair_coeff %s\n", $_);
}
printf(LAMMPS "\n");
}
- printf(LAMMPS "special_bonds charmm\n"); # invoke charmm
- printf(LAMMPS "thermo 1\n"); # set thermo style
+ 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 0.5\n\n"); # 0.5 ps time step
- printf(LAMMPS "minimize 0.0 0.0 1000 10000\n\n"); # take off the edge
+ 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
+ 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
+ if ($shake ne ""); # add water if present
printf(LAMMPS "velocity all create 0.0 12345678 dist uniform\n\n");
- printf(LAMMPS "restart 10 $project.restart1 $project.restart2\n");
- printf(LAMMPS "dump 1 all atom 10 $project.dump\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 "run 20\n"); # run for 20 time steps
+ printf(LAMMPS "thermo 100\n"); # set thermo style
+ printf(LAMMPS "run 1000\n"); # run for 1000 time steps
close(LAMMPS);
}