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