From c2e11dffa2eb23b253ebc6602bed58e7293431a2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 2 Oct 2016 20:33:20 -0400 Subject: [PATCH] import updated charmm2lammps.pl script from Rober Latour --- tools/ch2lmp/charmm2lammps.pl | 640 ++++++---------------------------- 1 file changed, 109 insertions(+), 531 deletions(-) diff --git a/tools/ch2lmp/charmm2lammps.pl b/tools/ch2lmp/charmm2lammps.pl index 0e39f589da..de8c29a646 100644 --- a/tools/ch2lmp/charmm2lammps.pl +++ b/tools/ch2lmp/charmm2lammps.pl @@ -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 and return the total number of $ident + + sub PSFGoto # goto $ident in { 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 @@ { # 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(" ", ); @@ -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 $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 = ; - 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); -