forked from lijiext/lammps
import updated charmm2lammps.pl script from Rober Latour
This commit is contained in:
Axel Kohlmeyer
parent
e2c7acabac
commit
c2e11dffa2
|
|
@ -28,43 +28,30 @@
|
|||
# 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
|
||||
# 20130508 Added 'CMAP crossterms' section at the end of the data file
|
||||
# 20131001 Added instructions in CMAP section to fix problem if 'ter'
|
||||
# is not designated in .pdb file to identify last amino acid
|
||||
#
|
||||
# 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
|
||||
#
|
||||
# General Many thanks to Paul S. Crozier for checking script validity
|
||||
# against his projects.
|
||||
#
|
||||
# ------------------------------------------------------------------------------
|
||||
# NOTE: This current version was modified by Xiaohu Hu (hux2@ornl.gov)
|
||||
# DATE: April, May 2009
|
||||
# Then finalized to complete addition of CMAP terms to data file
|
||||
# by Robert A. Latour, Clemson University (latourr@clemson.edu)
|
||||
# and Chris Lorenz, King's College (chris.lorenz@kcl.ac.uk)
|
||||
# DATE: May, 2013
|
||||
#
|
||||
# The original code of charmm2lammps.pl is modified
|
||||
#
|
||||
# 1. to fix the double-counting problem in 1-4 interaction associated
|
||||
# with the carbon-hydrate 6-rings containing systems. (See subroutine
|
||||
# DihedralCorrect6Ring, activated with the option "-cd")
|
||||
#
|
||||
# 2. to add a new section "CMAP" which is a list of the peptide
|
||||
# backbone dihedrals cross terms. A modifed LAMMPS version will be
|
||||
# needed to be able to use this feature. (See subroutine CharmmCmap,
|
||||
# activated with the option "-cmap")
|
||||
#
|
||||
# These subroutines are independent from the original charmm2lammps.pl, i.e.
|
||||
# the original code of charmm2lammps.pl is unchanged. The new routines only
|
||||
# evaluate and modify the output generated by the original charmm2lammps.pl.
|
||||
# -----------------------------------------------------------------------------
|
||||
|
||||
# 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(@_); # "@_" = argument passed to the subroutine
|
||||
my $name = shift(@_);
|
||||
|
||||
printf("Error: file %s not found\n", $name) if (!scalar(stat($name)));
|
||||
return !scalar(stat($name));
|
||||
|
|
@ -75,13 +62,11 @@
|
|||
{
|
||||
my $k = 0;
|
||||
my @dir = ("x", "y", "z");
|
||||
|
||||
# Modified by Xiaohu Hu, May 2009. Options "-cmap" and "-cdihedral" added
|
||||
my @options = ("-help", "-charmm", "-water", "-ions", "-center",
|
||||
my @options = ("-help", "-nohints", "-water", "-ions", "-center",
|
||||
"-quiet", "-pdb_ctrl", "-l", "-lx", "-ly", "-lz",
|
||||
"-border", "-ax", "-ay", "-az", "-cd", "-cmap");
|
||||
my @remarks = ("display this message",
|
||||
"add charmm types to LAMMPS data file",
|
||||
"-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",
|
||||
|
|
@ -95,19 +80,14 @@
|
|||
"rotation around x-axis",
|
||||
"rotation around y-axis",
|
||||
"rotation around z-axis",
|
||||
"use the 6-ring dihedral correction",
|
||||
"generate the CMAP section"
|
||||
);
|
||||
my $notes;
|
||||
|
||||
$program = "charmm2lammps";
|
||||
# $version = "1.8.2.5 beta"; # Modified by Xiaohu Hu, Dec 2009
|
||||
$version = "1.8.2.6 beta"; # Modified by Robert Latour & Chris Lorenz, May 2013
|
||||
# $year = "2009"; # Modified by Xiaohu Hu, April 2009
|
||||
$year = "2013"; # Modified by Robert Latour & Chris Lorenz, May 2013
|
||||
$cmap = 0; # Added by Xiaohu Hu, May 2009
|
||||
$cdihedral = 0; # Added by Xiaohu Hu, May 2009
|
||||
$add = 0;
|
||||
$version = "1.8.3";
|
||||
$year = "2016";
|
||||
$add = 1;
|
||||
$water_dens = 0;
|
||||
$ions = 0;
|
||||
$info = 1;
|
||||
|
|
@ -117,6 +97,7 @@
|
|||
$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";
|
||||
|
|
@ -141,7 +122,7 @@
|
|||
|
||||
foreach (@ARGV)
|
||||
{
|
||||
if (substr($_, 0, 1) eq "-") # if the first letter of the aguement is "-" = an option
|
||||
if (substr($_, 0, 1) eq "-")
|
||||
{
|
||||
my $k = 0;
|
||||
my @tmp = split("=");
|
||||
|
|
@ -152,24 +133,24 @@
|
|||
last if ($tmp[0] eq substr($_, 0 , length($tmp[0])));
|
||||
++$k;
|
||||
}
|
||||
$help = 1 if (!$k--);
|
||||
$add = 1 if (!$k--);
|
||||
$water_dens = ($tmp[1] ne "" ? $tmp[1] : 1) if (!$k--);
|
||||
$ion_molar = abs($tmp[1]) if (!$k);
|
||||
$ions = 1 if (!$k--);
|
||||
$center = 1 if (!$k--);
|
||||
$info = 0 if (!$k--);
|
||||
$pdb_ctrl = $switch if (!$k--);
|
||||
my $flag = $k--;
|
||||
$L[0] = abs($tmp[1]) if (!($flag && $k--));
|
||||
$L[1] = abs($tmp[1]) if (!($flag && $k--));
|
||||
$L[2] = abs($tmp[1]) if (!($flag && $k--));
|
||||
$border = abs($tmp[1]) if (!$k--);
|
||||
@R = M_Dot(M_Rotate(0, $tmp[1]), @R) if (!$k--);
|
||||
@R = M_Dot(M_Rotate(1, $tmp[1]), @R) if (!$k--);
|
||||
@R = M_Dot(M_Rotate(2, $tmp[1]), @R) if (!$k--);
|
||||
$cdihedral = 1 if (!$k--); # Added by Xiaohu Hu, May 2009
|
||||
$cmap = 1 if (!$k--); # Added by Xiaohu Hu, May 2009
|
||||
$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 = 1 if (!$k--);
|
||||
# print("Warning: ignoring unknown command line flag: $tmp[0]\n");
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
@ -180,7 +161,7 @@
|
|||
$water_dens = 1 if ($ions && !$water_dens);
|
||||
if (($k<2)||$help)
|
||||
{
|
||||
printf("\n%s v%s (c)%s by Pieter J. in \'t Veld for SNL\nwith 6-ring dihedral correction and CMAP added by X. Hu, 2009\n\n",
|
||||
printf("\n%s v%s (c)%s by Pieter J. in \'t Veld for SNL\n\n",
|
||||
$program, $version, $year);
|
||||
printf("Usage:\n %s.pl [-option[=#] ..] forcefield project\n\n",$program);
|
||||
printf("Options:\n");
|
||||
|
|
@ -191,7 +172,7 @@
|
|||
printf("\nNotes:\n%s\n", $notes);
|
||||
exit(-1);
|
||||
}
|
||||
else { printf("%s v%s (c)%s\n\n", $program, $version, $year) if ($info); }
|
||||
else { printf("\n%s v%s (c)%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")
|
||||
|
|
@ -330,7 +311,7 @@
|
|||
|
||||
sub PSFConnectivity
|
||||
{
|
||||
my $n = PSFGoto(bonds); # $n = the total number of bonds
|
||||
my $n = PSFGoto(bonds);
|
||||
|
||||
return if (scalar(@nconnect));
|
||||
printf("Info: creating connectivity\n") if ($info);
|
||||
|
|
@ -342,6 +323,7 @@
|
|||
}
|
||||
}
|
||||
|
||||
|
||||
sub PSFDihedrals # hack to accomodate
|
||||
{ # LAMMPS' way of calc
|
||||
$idihedral = 0; # LJ 1-4 interactions
|
||||
|
|
@ -372,6 +354,7 @@
|
|||
return $ndihedral;
|
||||
}
|
||||
|
||||
|
||||
sub CreatePSFIndex # make an index of id
|
||||
{ # locations
|
||||
my @psf_ids = ("!NATOM","!NBOND:","!NTHETA:","!NPHI:","!NIMPHI:");
|
||||
|
|
@ -391,24 +374,26 @@
|
|||
}
|
||||
}
|
||||
|
||||
sub PSFGoto # goto $ident in <PSF> and return the total number of $ident
|
||||
|
||||
sub PSFGoto # goto $ident in <PSF>
|
||||
{
|
||||
CreatePSFIndex() if (!scalar(%PSFIndex));
|
||||
my $id = shift(@_);
|
||||
my @n = split(" ", $PSFIndex{$id}); # = 1 if the $ident is found
|
||||
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); # set file-handle position to the EOF
|
||||
seek(PSF, 0, SEEK_END);
|
||||
return -1;
|
||||
}
|
||||
seek(PSF, $n[0], SEEK_SET);
|
||||
return $n[1];
|
||||
}
|
||||
|
||||
|
||||
sub PSFGet
|
||||
{
|
||||
if ($dihedral_flag)
|
||||
|
|
@ -426,7 +411,7 @@
|
|||
}
|
||||
|
||||
|
||||
sub PSFWrite
|
||||
sub PSFWrite
|
||||
{
|
||||
my $items = shift(@_);
|
||||
my $n = $items;
|
||||
|
|
@ -542,11 +527,14 @@
|
|||
|
||||
sub Markers
|
||||
{
|
||||
my %markers;
|
||||
my $n = 0;
|
||||
|
||||
foreach ("NONBONDED", "BONDS", "ANGLES", "DIHEDRALS", "IMPROPER") {
|
||||
$markers{$_} = $n++; }
|
||||
my %markers = (
|
||||
NONBONDED => '0',
|
||||
BONDS => '1',
|
||||
ANGLES => '2',
|
||||
DIHEDRALS => '3',
|
||||
IMPROPERS => '4',
|
||||
IMPROPER => '4'
|
||||
);
|
||||
return %markers;
|
||||
}
|
||||
|
||||
|
|
@ -636,7 +624,7 @@
|
|||
{ # <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;
|
||||
|
|
@ -711,7 +699,7 @@
|
|||
|
||||
sub SetScreeningFactor # set screening factor
|
||||
{
|
||||
my $id = shift(@_); # @_ has the form: ($variable, float no.)
|
||||
my $id = shift(@_);
|
||||
my $value = shift(@_);
|
||||
my $new = "";
|
||||
|
||||
|
|
@ -725,6 +713,7 @@
|
|||
$parms[$id] = $new;
|
||||
}
|
||||
|
||||
|
||||
sub CorrectDihedralParameters
|
||||
{
|
||||
my $n = PSFGoto(dihedrals);
|
||||
|
|
@ -735,20 +724,15 @@
|
|||
my $first;
|
||||
my $last;
|
||||
|
||||
for (my $i=0; $i<$n; ++$i) # loop over all dihedrals
|
||||
for (my $i=0; $i<$n; ++$i)
|
||||
{
|
||||
my @bonded = PSFGet(4);
|
||||
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 the condition "$id2 = $hash{$hash_id = $first." ".$last}" is false, which means
|
||||
# the id isn't in the hash, then it will be added, otherwise means this condition is
|
||||
# true and the
|
||||
if (($id2 = $hash{$hash_id = $first." ".$last}) eq "")
|
||||
{
|
||||
$hash{$hash_id} = $id1; # add id to hash
|
||||
|
|
@ -1037,7 +1021,7 @@
|
|||
|
||||
sub WriteLAMMPSHeader # print lammps header
|
||||
{
|
||||
printf(LAMMPS "Created by $program v$version on %s\n", `date`);
|
||||
printf(LAMMPS "LAMMPS data file. CGCMM style. atom_style full. Created by $program v$version on %s\n", `date`);
|
||||
printf(LAMMPS "%12d atoms\n", $natoms);
|
||||
printf(LAMMPS "%12d bonds\n", $nbonds);
|
||||
printf(LAMMPS "%12d angles\n", $nangles);
|
||||
|
|
@ -1191,7 +1175,7 @@
|
|||
|
||||
CRDGoto(atoms);
|
||||
$net_charge = 0;
|
||||
printf(LAMMPS "Atoms\n\n") if ($n>0);
|
||||
printf(LAMMPS "Atoms # full\n\n") if ($n>0);
|
||||
for (my $i=0; $i<$n; ++$i)
|
||||
{
|
||||
my @crd = $pdb ? NextPDB2CRD() : split(" ", <CRD>);
|
||||
|
|
@ -1225,6 +1209,7 @@
|
|||
{
|
||||
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;
|
||||
|
||||
|
|
@ -1238,7 +1223,7 @@
|
|||
@parms = Delete(1, @parms) if ($mode==3);
|
||||
}
|
||||
return 0 if (!scalar(@parms));
|
||||
printf(LAMMPS "%s Coeffs\n\n", $header);
|
||||
printf(LAMMPS "%s Coeffs # %s\n\n", $header, $hint);
|
||||
for (my $i=0; $i<scalar(@parms); ++$i)
|
||||
{
|
||||
if ($parms[$i] ne "")
|
||||
|
|
@ -1418,6 +1403,11 @@
|
|||
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);
|
||||
|
|
@ -1449,27 +1439,22 @@
|
|||
printf(LAMMPS "# Created by $program v$version on %s\n", `date`);
|
||||
printf(LAMMPS "units real\n"); # general
|
||||
printf(LAMMPS "neigh_modify delay 2 every 1\n\n");
|
||||
|
||||
printf(LAMMPS "boundary p p p\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 "# if cutoffs to be used for electrostatics, use pair_style lj/charmmfsw/coul/charmmfsh\n");
|
||||
printf(LAMMPS "# and delete or comment out kspace_style\n");
|
||||
printf(LAMMPS "pair_style lj/charmm/coul/long 8 12\n");
|
||||
printf(LAMMPS "pair_modify mix arithmetic\n");
|
||||
printf(LAMMPS "kspace_style pppm 1e-4\n\n");
|
||||
|
||||
|
||||
if ($cmap) {
|
||||
printf(LAMMPS "# Modify following lines to provide directory path cmap.data and 'project'.data files\n");
|
||||
printf(LAMMPS "fix cmap all cmap /directoryPath/.../cmap36.data\n");
|
||||
printf(LAMMPS "fix_modify cmap energy yes\n");
|
||||
printf(LAMMPS "read_data /directoryPath/.../$project.data fix cmap crossterm CMAP\n\n");
|
||||
printf(LAMMPS "read_data $project.data fix cmap crossterm CMAP\n\n");
|
||||
}else{
|
||||
printf(LAMMPS "# Modify following line to provide directory path for 'project'.data file\n");
|
||||
printf(LAMMPS "read_data /directoryPath/.../$project.data\n\n"); # read data
|
||||
printf(LAMMPS "read_data $project.data\n\n"); # read data
|
||||
}
|
||||
|
||||
if ($coefficients ne "") # corrected coeffs
|
||||
|
|
@ -1481,15 +1466,16 @@
|
|||
printf(LAMMPS "\n");
|
||||
}
|
||||
printf(LAMMPS "special_bonds charmm\n"); # invoke charmm
|
||||
printf(LAMMPS "thermo 1\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 "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 "thermo 1\n"); # set thermo style
|
||||
printf(LAMMPS "thermo_style multi\n");
|
||||
printf(LAMMPS "timestep 0.5\n\n"); # 0.5 ps time step
|
||||
printf(LAMMPS "restart 10 $project.restart1 $project.restart2\n");
|
||||
printf(LAMMPS "dump 1 all atom 10 $project.dump\n");
|
||||
printf(LAMMPS "dump_modify 1 image yes scale yes\n\n");
|
||||
|
|
@ -1497,435 +1483,32 @@
|
|||
close(LAMMPS);
|
||||
}
|
||||
|
||||
# ------------------ DESCRIPTION: sub DihedralCorrect6Ring ------------------ #
|
||||
# #
|
||||
# This subroutine is a subsequent correction of the dihedral 1-4 non-bonded #
|
||||
# interaction scaling factor in the LAMMPS data file. The problem occurs in #
|
||||
# dealing with carbohydrate systems, when some dihedrals outside of a 6-ring #
|
||||
# are incorrectly assigned to the same dihedral type as dihedrals within a #
|
||||
# 6-ring. Thus, those dihedrals outside of a 6-ring will be treated with a #
|
||||
# 1-4 non-bonded interaction scaling factor of 0.5 instead of 1. #
|
||||
# #
|
||||
# By: Xiaohu Hu (hux2@ornl.gov) #
|
||||
# #
|
||||
# April 2009 #
|
||||
# #
|
||||
# --------------------------------------------------------------------------- #
|
||||
# ----------------------- 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 DihedralCorrect6Ring
|
||||
{
|
||||
print "\nINITIATE DIHEDRAL CORRECTION...\n\n";
|
||||
## Variables for general data
|
||||
# arrays
|
||||
my @temp_data1;
|
||||
my @temp_data2;
|
||||
my @dihedral_coeffs;
|
||||
my @dihedral_list;
|
||||
my @raw_data;
|
||||
my @temp1;
|
||||
my @temp2;
|
||||
my @temp3;
|
||||
my @ring_dihe_list;
|
||||
my @ring_dihetype_list;
|
||||
|
||||
# integers
|
||||
my $counter1 = 0;
|
||||
my $counter2 = 0;
|
||||
my $net_ndihedral = 0;
|
||||
my $net_ndihedral_types = 0;
|
||||
my $splice_onset_dihe = 0;
|
||||
my $splice_onset_dihe_coeff = 0;
|
||||
my $n;
|
||||
my $ntyp;
|
||||
|
||||
# matrix
|
||||
my @dihedral_matrix;
|
||||
my @dihedral_coeff_matrix;
|
||||
my @new_dihedral_coeff_matrix;
|
||||
my @new_dihedral_matrix;
|
||||
|
||||
## Variables for dihedral matrix data
|
||||
# "Dihedral Coeffs"
|
||||
my $dihedral_type_c;
|
||||
my $amp;
|
||||
my $period;
|
||||
my $equ_val;
|
||||
my $scaling;
|
||||
|
||||
# "Dihedrals"
|
||||
my $dihedral_no;
|
||||
my $dihedral_type_l;
|
||||
my $atom_1_no;
|
||||
my $atom_2_no;
|
||||
my $atom_3_no;
|
||||
my $atom_4_no;
|
||||
|
||||
# ------------- Reread the previously generated LAMMPS data file --------------------
|
||||
|
||||
open(LAMMPS, "< $project.data") or
|
||||
die "\"sub DihedralCorrect6Ring\" cannot open \"$project.data!\n";
|
||||
print "Analysing \"$project.data\"...\n\n";
|
||||
@raw_data = <LAMMPS>;
|
||||
close(LAMMPS);
|
||||
|
||||
# Find the number of dihedrals and dihedral types and the sections "Dihedrals"
|
||||
# and "Dihedral Coeffs" in LAMMPS data file
|
||||
foreach $line (@raw_data) {
|
||||
$counter1++;
|
||||
chomp($line);
|
||||
if ($line =~ m/dihedrals/) {
|
||||
($n,$string) = split(" ",$line);
|
||||
print "$n dihedrals\n";
|
||||
}
|
||||
if ($line =~ m/dihedral types/) {
|
||||
($ntyp,$string) = split(" ",$line);
|
||||
print "$ntyp dihedral types\n";
|
||||
}
|
||||
if ($line =~ m/Dihedral Coeffs/) {
|
||||
print "Section \"Dihedrals\" found\n";
|
||||
$splice_onset_dihe_coeff = $counter1;
|
||||
}
|
||||
if ($line =~ m/Dihedrals/) {
|
||||
print "Section \"Dihedral Coeffs\" found\n";
|
||||
$splice_onset_dihe = $counter1;
|
||||
}
|
||||
}
|
||||
if ($splice_onset_dihe_coeff == 0 and $splice_onset_dihe == 0) {
|
||||
print "\nNo dihedral data in \"$project.data\", no dihedral correction necessary\n";
|
||||
return;
|
||||
}
|
||||
elsif ($splice_onset_dihe_coeff == 0 or $splice_onset_dihe == 0) {
|
||||
print "\nDihedral data incomplete! Dihedral correction impossible\n";
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
# --------------- Transform the raw data into matrices -----------------------
|
||||
|
||||
@temp_data1 = @temp_data2 = @raw_data;
|
||||
|
||||
#Store the section "Dihedral Coeffs" into a new list
|
||||
@dihedral_coeffs = splice(@temp_data1,$splice_onset_dihe_coeff,$ntyp+1);
|
||||
|
||||
# Transfer the data from "@dihedral_coeffs" (an array of strings)
|
||||
# into "@dihedral_coeff_matrix" (a 6 x $n matrix of integers)
|
||||
for (@dihedral_coeffs) {
|
||||
($dihedral_type_c,
|
||||
$amp,
|
||||
$period,
|
||||
$equ_val,
|
||||
$scaling) = split(" ");
|
||||
|
||||
push(@dihedral_coeff_matrix,
|
||||
[$dihedral_type_c,
|
||||
$amp,
|
||||
$period,
|
||||
$equ_val,
|
||||
$scaling]);
|
||||
}
|
||||
|
||||
# Store the section "Dihedrals" into a new list
|
||||
@dihedral_list = splice(@temp_data2,$splice_onset_dihe,$n+1);
|
||||
|
||||
# Transfer the data from "@dihedral_list" (an array of strings)
|
||||
# into "@dihedral_matrix" (a 6 x $n matrix of integers)
|
||||
for (@dihedral_list) {
|
||||
($dihedral_no,
|
||||
$dihedral_type_l,
|
||||
$atom_no_1,
|
||||
$atom_no_2,
|
||||
$atom_no_3,
|
||||
$atom_no_4) = split(" ");
|
||||
|
||||
push(@dihedral_matrix,
|
||||
[$dihedral_no,
|
||||
$dihedral_type_l,
|
||||
$atom_no_1,
|
||||
$atom_no_2,
|
||||
$atom_no_3,
|
||||
$atom_no_4]);
|
||||
}
|
||||
|
||||
for (my $i = 1; $i < $n; $i++) {
|
||||
my $cur_type = ${$dihedral_matrix[$i]}[1];
|
||||
if ($cur_type == 1 or
|
||||
$cur_type == 16 or
|
||||
$cur_type == 34 or
|
||||
$cur_type == 46 or
|
||||
$cur_type == 58 or
|
||||
$cur_type == 64) {
|
||||
push(@list,$cur_type);
|
||||
}
|
||||
}
|
||||
@list = sort {$a <=> $b} @list;
|
||||
#print "@list\n";
|
||||
#print "Total: $#list\n";
|
||||
|
||||
# ------------------ Reformat the matrices -------------------------------
|
||||
|
||||
# Loop the dihedral coefficient matrix and throw away all elements
|
||||
# with a zero scaling factor (= entries with periodicity != 1)
|
||||
for (my $i = 1; $i < $ntyp; $i++) {
|
||||
my $current_sf = ${$dihedral_coeff_matrix[$i]}[4];
|
||||
if ($current_sf != 0) {
|
||||
push(@new_dihedral_coeff_matrix,
|
||||
[${$dihedral_coeff_matrix[$i]}[0],
|
||||
${$dihedral_coeff_matrix[$i]}[1],
|
||||
${$dihedral_coeff_matrix[$i]}[2],
|
||||
${$dihedral_coeff_matrix[$i]}[3],
|
||||
${$dihedral_coeff_matrix[$i]}[4]]);
|
||||
$net_ndihedral_types++;
|
||||
}
|
||||
}
|
||||
|
||||
# Remove the duplicated dihedrals from the @dihedral_matrix and save
|
||||
# results into @new_dihedral_matrix
|
||||
push(@new_dihedral_matrix,
|
||||
[${$dihedral_matrix[1]}[0],
|
||||
${$dihedral_matrix[1]}[1],
|
||||
${$dihedral_matrix[1]}[2],
|
||||
${$dihedral_matrix[1]}[3],
|
||||
${$dihedral_matrix[1]}[4],
|
||||
${$dihedral_matrix[1]}[5]]);
|
||||
$net_ndihedral = 1;
|
||||
for (my $i = 2; $i < $n; $i++) {
|
||||
if (${$dihedral_matrix[$i]}[2] != ${$dihedral_matrix[$i-1]}[2] or
|
||||
${$dihedral_matrix[$i]}[3] != ${$dihedral_matrix[$i-1]}[3] or
|
||||
${$dihedral_matrix[$i]}[4] != ${$dihedral_matrix[$i-1]}[4] or
|
||||
${$dihedral_matrix[$i]}[5] != ${$dihedral_matrix[$i-1]}[5]) {
|
||||
push(@new_dihedral_matrix,
|
||||
[${$dihedral_matrix[$i]}[0],
|
||||
${$dihedral_matrix[$i]}[1],
|
||||
${$dihedral_matrix[$i]}[2],
|
||||
${$dihedral_matrix[$i]}[3],
|
||||
${$dihedral_matrix[$i]}[4],
|
||||
${$dihedral_matrix[$i]}[5]]);
|
||||
$net_ndihedral++;
|
||||
}
|
||||
}
|
||||
|
||||
# Print out the matrix
|
||||
# my $ref_line;
|
||||
# my $column;
|
||||
# print "New dihedral list:\n";
|
||||
# foreach $ref_line (@new_dihedral_matrix) {
|
||||
# foreach $column (@$ref_line) {
|
||||
# print "$column "
|
||||
# }
|
||||
# print "\n";
|
||||
# }
|
||||
|
||||
# --------------- Seach for the wrong scaling factors ----------------------------
|
||||
|
||||
# Loop through the dihedral matrix to determine the dihedrals within
|
||||
# a 6-ring. Note there is some duplication in this approach due to
|
||||
# processing already processed dihedrals. This will be taken care of
|
||||
# later.
|
||||
#
|
||||
# NOTE: @ringlist contains the dihedrals types which corresponds to
|
||||
# dihedrals within a 6-ring system.
|
||||
#
|
||||
my $n6ring = 0;
|
||||
for (my $i = 1; $i < $net_ndihedral; $i++) {
|
||||
my $current_i = ${$new_dihedral_matrix[$i]}[2];
|
||||
my $current_l = ${$new_dihedral_matrix[$i]}[5];
|
||||
for (my $j = $i + 1; $j < $net_ndihedral; $j++) {
|
||||
if ( ($current_i == ${$new_dihedral_matrix[$j]}[2] and
|
||||
$current_l == ${$new_dihedral_matrix[$j]}[5]) or
|
||||
($current_i == ${$new_dihedral_matrix[$j]}[5] and
|
||||
$current_l == ${$new_dihedral_matrix[$j]}[2]) ) {
|
||||
push(@temp2,${$new_dihedral_matrix[$i]}[0]);
|
||||
push(@temp2,${$new_dihedral_matrix[$j]}[0]);
|
||||
push(@temp3,${$new_dihedral_matrix[$i]}[1]);
|
||||
push(@temp3,${$new_dihedral_matrix[$j]}[1]);
|
||||
$n6ring++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if ($n6ring == 0) {
|
||||
print "\nNo dihedrals within 6-ring structure found. No correction necessary.\n";
|
||||
return;
|
||||
}
|
||||
|
||||
# Sort the lists according to the numerical order
|
||||
@ring_dihe_list = sort {$a <=> $b} @temp2;
|
||||
@ring_dihetype_list = sort {$a <=> $b} @temp3;
|
||||
|
||||
# Remove the dups
|
||||
my %seen1;
|
||||
for (my $i = 0; $i <= $#ring_dihe_list;) {
|
||||
splice @ring_dihe_list, --$i, 1
|
||||
if $seen1 {$ring_dihe_list[$i++]}++;
|
||||
}
|
||||
|
||||
my %seen2;
|
||||
for (my $i = 0; $i <= $#ring_dihetype_list;) {
|
||||
splice @ring_dihetype_list, --$i, 1
|
||||
if $seen1 {$ring_dihetype_list[$i++]}++;
|
||||
}
|
||||
|
||||
print "6-ring dihedral types:\n";
|
||||
print "@ring_dihetype_list\n";
|
||||
print "6-ring dihedrals:\n";
|
||||
print "@ring_dihe_list\n";
|
||||
|
||||
# Locate the wrong dihedrals in the dihedral list. Criteria to be wrong:
|
||||
# if the type of the dihedral is equal to one from the list @ring_dihetype_list
|
||||
# but not one from the @ring_dihe_list
|
||||
my @errorlist;
|
||||
my @errortypes;
|
||||
my @raw_errortypes;
|
||||
my $dihe_flag = 0;
|
||||
|
||||
for (my $i = 1; $i < $n; $i++) {
|
||||
my $cur_dihe = ${dihedral_matrix[$i]}[0];
|
||||
my $cur_type = ${dihedral_matrix[$i]}[1];
|
||||
|
||||
for (my $j = 0; $j <= $#ring_dihetype_list; $j++) {
|
||||
if ($cur_type == $ring_dihetype_list[$j]) {
|
||||
for (my $k = 0; $k <= $#ring_dihe_list; $k++) {
|
||||
if ($cur_dihe == $ring_dihe_list[$k]) {
|
||||
$dihe_flag++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if ($dihe_flag == 0) {
|
||||
for (my $j = 0; $j <= $#ring_dihetype_list; $j++) {
|
||||
if ($cur_type == $ring_dihetype_list[$j]) {
|
||||
push(@errorlist,$cur_dihe);
|
||||
push(@raw_errortypes,$cur_type);
|
||||
$counter2++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
$dihe_flag = 0;
|
||||
}
|
||||
|
||||
if ($counter2 == 0) {
|
||||
print "\nNo mis-assigned dihedrals found. No correction necessary.\n";
|
||||
return;
|
||||
}
|
||||
|
||||
print "mis-assigned dihedral/s found: $counter2\n";
|
||||
print "@errorlist\n";
|
||||
# print "@raw_errortypes\n";
|
||||
|
||||
# Sort the @errortypes and remove dups
|
||||
@errortypes = sort {$a <=> $b} @raw_errortypes;
|
||||
my %seen4;
|
||||
for (my $i = 0; $i <= $#errortypes;) {
|
||||
splice @errortypes, --$i, 1
|
||||
if $seen4 {$errortypes[$i++]}++;
|
||||
}
|
||||
|
||||
print "associated dihedral type/s\n";
|
||||
print "@errortypes\n";
|
||||
|
||||
# ------------ Add new dihedral types for the mis-assigned dihedrals -----------
|
||||
|
||||
print "\nWriting \"corrected_$project.data\"...\n\n";
|
||||
open(REWRITE,"> corrected_$project.data")
|
||||
or die "Can not write file \"corrected_$project.data\"!\n";
|
||||
|
||||
my $counter3 = 0;
|
||||
my $fix_start = $splice_onset_dihe_coeff + $ntyp + 1;
|
||||
my $new_ntyp = $ntyp + $#errortypes +1;
|
||||
|
||||
foreach $line (@raw_data) {
|
||||
# update header information
|
||||
if ($line =~ m/dihedral types/) {
|
||||
$line =~ s/$ntyp/$new_ntyp/;
|
||||
}
|
||||
|
||||
$counter3++;
|
||||
printf(REWRITE "$line\n");
|
||||
if ($counter3 == $fix_start) { last }
|
||||
}
|
||||
|
||||
my @newtypes;
|
||||
for (my $i = 0; $i <= $#errortypes; $i++) {
|
||||
$ntyp++;
|
||||
printf(REWRITE "%8d %10.7g %10d %10d 1\n",
|
||||
$ntyp,
|
||||
${$dihedral_coeff_matrix[$errortypes[$i]]}[1],
|
||||
${$dihedral_coeff_matrix[$errortypes[$i]]}[2],
|
||||
${$dihedral_coeff_matrix[$errortypes[$i]]}[3]);
|
||||
print "New dihedral type $ntyp added\n";
|
||||
push(@newtypes,$ntyp);
|
||||
}
|
||||
|
||||
# ------------ Assign the wrong dihedrals with the new types ----------------
|
||||
|
||||
print "\nCorrecting the mis-assigned dihedrals...\n\n";
|
||||
|
||||
printf(REWRITE "$raw_data[$splice_onset_dihe - 2]\n");
|
||||
printf(REWRITE "$raw_data[$splice_onset_dihe - 1]\n");
|
||||
|
||||
# Overwrite the wrong dihedrals in the dihedral list (section "Dihedrals")
|
||||
my $counter3 = -1;
|
||||
my $write_flag = 0;
|
||||
my @temp_line;
|
||||
|
||||
for ($j = $splice_onset_dihe; $j <= $#raw_data; $j++) {
|
||||
$counter3++;
|
||||
for (my $i = 0; $i <= $#errorlist; $i++) {
|
||||
if ($counter3 == $errorlist[$i]) {
|
||||
@temp_line = split(" ",$raw_data[$j]);
|
||||
for (my $k = 0; $k <= $#errortypes; $k++) {
|
||||
if ($temp_line[1] == $errortypes[$k]) {
|
||||
printf(REWRITE "%8d %7d %7d %7d %7d %7d\n",
|
||||
${$dihedral_matrix[$errorlist[$i]]}[0],
|
||||
$newtypes[$k],
|
||||
${$dihedral_matrix[$errorlist[$i]]}[2],
|
||||
${$dihedral_matrix[$errorlist[$i]]}[3],
|
||||
${$dihedral_matrix[$errorlist[$i]]}[4],
|
||||
${$dihedral_matrix[$errorlist[$i]]}[5]);
|
||||
print "Dihedral No. ${$dihedral_matrix[$errorlist[$i]]}[0] corrected\n";
|
||||
$write_flag = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if ($write_flag == 0) { printf(REWRITE "$raw_data[$j]\n");}
|
||||
$write_flag = 0;
|
||||
}
|
||||
|
||||
print "\nDONE!\n";
|
||||
|
||||
close(REWRITE);
|
||||
|
||||
# End of the subroutine
|
||||
}
|
||||
|
||||
# ----------------------- 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) #
|
||||
# #
|
||||
# By: Xiaohu Hu (hux2@ornl.gov) #
|
||||
# May 2009 #
|
||||
# #
|
||||
# Modified May 2013 by: Robert Latour (latourr@clemson.edu) 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
|
||||
{
|
||||
sub CharmmCmap
|
||||
{
|
||||
print "\nINITIATING CHARMM CMAP SUBROUTINE...\n\n";
|
||||
|
||||
# Reread and analyse $project.data
|
||||
|
|
@ -2458,10 +2041,7 @@ sub CharmmCmap
|
|||
printf(REWRITE "\n");
|
||||
}
|
||||
close(REWRITE);
|
||||
|
||||
print "\nDone!\n\n";
|
||||
|
||||
# End of the CharmmCmap subroutine
|
||||
}
|
||||
|
||||
# main
|
||||
|
|
@ -2470,6 +2050,4 @@ sub CharmmCmap
|
|||
WriteData();
|
||||
WriteLAMMPSInput();
|
||||
printf("Info: conversion complete\n\n") if ($info);
|
||||
DihedralCorrect6Ring() if ($cdihedral);
|
||||
CharmmCmap() if ($cmap);
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue