diff --git a/doc/src/Manual.txt b/doc/src/Manual.txt index 441e9cc8f0..cd8750bd84 100644 --- a/doc/src/Manual.txt +++ b/doc/src/Manual.txt @@ -1,7 +1,7 @@ LAMMPS Users Manual - + @@ -21,7 +21,7 @@

LAMMPS Documentation :c,h3 -21 Sep 2016 version :c,h4 +26 Sep 2016 version :c,h4 Version info: :h4 diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt index 3e18037a83..91ec679360 100644 --- a/doc/src/Section_commands.txt +++ b/doc/src/Section_commands.txt @@ -501,6 +501,7 @@ USER-INTEL, k = KOKKOS, o = USER-OMP, t = OPT. "bond/create"_fix_bond_create.html, "bond/swap"_fix_bond_swap.html, "box/relax"_fix_box_relax.html, +"cmap"_fix_cmap.html, "controller"_fix_controller.html, "deform (k)"_fix_deform.html, "deposit"_fix_deposit.html, diff --git a/doc/src/Section_tools.txt b/doc/src/Section_tools.txt index 3527d228de..983be30f19 100644 --- a/doc/src/Section_tools.txt +++ b/doc/src/Section_tools.txt @@ -107,9 +107,10 @@ The ch2lmp sub-directory contains tools for converting files back-and-forth between the CHARMM MD code and LAMMPS. They are intended to make it easy to use CHARMM as a builder and as a -post-processor for LAMMPS. Using charmm2lammps.pl, you can convert an -ensemble built in CHARMM into its LAMMPS equivalent. Using -lammps2pdb.pl you can convert LAMMPS atom dumps into pdb files. +post-processor for LAMMPS. Using charmm2lammps.pl, you can convert a +PDB file with associated CHARMM info, including CHARMM force field +data, into its LAMMPS equivalent. Using lammps2pdb.pl you can convert +LAMMPS atom dumps into PDB files. See the README file in the ch2lmp sub-directory for more information. diff --git a/doc/src/fix_cmap.txt b/doc/src/fix_cmap.txt new file mode 100644 index 0000000000..924e6f5022 --- /dev/null +++ b/doc/src/fix_cmap.txt @@ -0,0 +1,132 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +fix cmap command :h3 + +[Syntax:] + +fix ID group-ID cmap filename :pre + +ID, group-ID are documented in "fix"_fix.html command +cmap = style name of this fix command +filename = force-field file with CMAP coefficients :ul + +[Examples:] + +fix myCMAP all cmap ../potentials/cmap36.data +read_data proteinX.data fix myCMAP crossterm CMAP +fix_modify myCMAP energy yes :pre + +[Description:] + +This command enables CMAP crossterms to be added to simulations which +use the CHARMM force field. These are relevant for any CHARMM model +of a peptide or protein sequences that is 3 or more amino-acid +residues long; see "(Buck)"_#Buck and "(Brooks)"_#Brooks for details, +including the analytic energy expressions for CMAP interactions. The +CMAP crossterms add additional potential energy contributions to pairs +of overlapping phi-psi dihedrals of amino-acids, which are important +to properly represent their conformational behavior. + +The examples/cmap directory has a sample input script and data file +for a small peptide, that illustrates use of the fix cmap command. + +As in the example above, this fix should be used before reading a data +file that contains a listing of CMAP interactions. The {filename} +specified should contain the CMAP parameters for a particular version +of the CHARMM force field. Two such files are including in the +lammps/potentials directory: charmm22.cmap and charmm36.cmap. + +The data file read by the "read_data" must contain the topology of all +the CMAP interactions, similar to the topology data for bonds, angles, +dihedrals, etc. Specically it should have a line like this +in its header section: + +N crossterms :pre + +where N is the number of CMAP crossterms. It should also have a section +in the body of the data file like this with N lines: + +CMAP :pre + + 1 1 8 10 12 18 20 + 2 5 18 20 22 25 27 + ... + N 3 314 315 317 318 330 :pre + +The first column is an index from 1 to N to enumerate the CMAP terms; +it is ignored by LAMMPS. The 2nd column is the "type" of the +interaction; it is an index into the CMAP force field file. The +remaining 5 columns are the atom IDs of the atoms in the two 4-atom +dihedrals that overlap to create the CMAP 5-body interaction. Note +that the "crossterm" and "CMAP" keywords for the header and body +sections match those specified in the read_data command following the +data file name; see the "read_data"_doc/read_data.html doc page for +more details. + +A data file containing CMAP crossterms can be generated from a PDB +file using the charmm2lammps.pl script in the tools/ch2lmp directory +of the LAMMPS distribution. The script must be invoked with the +optional "-cmap" flag to do this; see the tools/ch2lmp/README file for +more information. + +The potential energy associated with CMAP interactions can be output +as described below. It can also be included in the total potential +energy of the system, as output by the +"thermo_style"_thermo_style.html command, if the "fix_modify +energy"_fix_modify.html command is used, as in the example above. See +the note below about how to include the CMAP energy when performing an +"energy minimization"_minimize.html. + +:line + +[Restart, fix_modify, output, run start/stop, minimize info:] + +No information about this fix is written to "binary restart +files"_restart.html. + +The "fix_modify"_fix_modify.html {energy} option is supported by this +fix to add the potential "energy" of the CMAP interactions system's +potential energy as part of "thermodynamic output"_thermo_style.html. + +This fix computes a global scalar which can be accessed by various +"output commands"_Section_howto.html#howto_15. The scalar is the +potential energy discussed above. The scalar value calculated by this +fix is "extensive". + +No parameter of this fix can be used with the {start/stop} keywords of +the "run"_run.html command. + +The forces due to this fix are imposed during an energy minimization, +invoked by the "minimize"_minimize.html command. + +NOTE: If you want the potential energy associated with the CMAP terms +forces to be included in the total potential energy of the system (the +quantity being minimized), you MUST enable the +"fix_modify"_fix_modify.html {energy} option for this fix. + +[Restrictions:] + +This fix can only be used if LAMMPS was built with the MOLECULE +package (which it is by default). See the "Making +LAMMPS"_Section_start.html#start_3 section for more info on packages. + +[Related commands:] + +"fix_modify"_fix_modify.html, "read_data"_read_data.html + +[Default:] none + +:line + +(Buck) +Buck, Bouguet-Bonnet, Pastor, MacKerell Jr., Biophys J, 90, L36 +(2006). + +(Brooks) +Brooks, Brooks, MacKerell Jr., J Comput Chem, 30, 1545 (2009). diff --git a/examples/README b/examples/README index 400d69aa2d..462fa93e12 100644 --- a/examples/README +++ b/examples/README @@ -61,6 +61,7 @@ sub-directories: accelerate: use of all the various accelerator packages balance: dynamic load balancing, 2d system body: body particles, 2d system +cmap: CMAP 5-body contributions to CHARMM force field colloid: big colloid particles in a small particle solvent, 2d system coreshell: adiabatic core/shell model comb: models using the COMB potential diff --git a/examples/cmap/charmm22.cmap b/examples/cmap/charmm22.cmap new file mode 100644 index 0000000000..782dcc5dbf --- /dev/null +++ b/examples/cmap/charmm22.cmap @@ -0,0 +1,1022 @@ +# Title: charmm correction map + +# alanine map, type 1 + +# -180 +0.126790 0.768700 0.971260 1.250970 2.121010 +2.695430 2.064440 1.764790 0.755870 -0.713470 +0.976130 -2.475520 -5.455650 -5.096450 -5.305850 +-3.975630 -3.088580 -2.784200 -2.677120 -2.646060 +-2.335350 -2.010440 -1.608040 -0.482250 + +# -165 +-0.802290 1.377090 1.577020 1.872290 2.398990 +2.461630 2.333840 1.904070 1.061460 0.518400 +-0.116320 -3.575440 -5.284480 -5.160310 -4.196010 +-3.276210 -2.715340 -1.806200 -1.101780 -1.210320 +-1.008810 -0.637100 -1.603360 -1.776870 + +# -150 +-0.634810 1.156210 1.624350 2.047200 2.653910 +2.691410 2.296420 1.960450 1.324930 2.038290 +-1.151510 -3.148610 -4.058280 -4.531850 -3.796370 +-2.572090 -1.727250 -0.961410 -0.282910 -0.479120 +-1.039340 -1.618060 -1.725460 -1.376360 + +# -135 +0.214000 1.521370 1.977440 2.377950 2.929470 +2.893410 2.435810 2.162970 1.761500 1.190090 +-1.218610 -2.108900 -2.976100 -3.405340 -2.768440 +-1.836030 -0.957950 0.021790 -0.032760 -0.665880 +-1.321170 -1.212320 -0.893170 -0.897040 + +# -120 +0.873950 1.959160 2.508990 2.841100 3.698960 +3.309330 2.614300 2.481720 2.694660 1.082440 +-0.398320 -1.761800 -2.945110 -3.294690 -2.308300 +-0.855480 -0.087320 0.439040 0.691880 -0.586330 +-1.027210 -0.976640 -0.467580 0.104020 + +# -105 +1.767380 2.286650 2.818030 3.065500 3.370620 +3.397440 2.730310 2.878790 2.542010 1.545240 +-0.092150 -1.694440 -2.812310 -2.802430 -1.856360 +-0.306240 -0.122440 0.444680 0.810150 -0.058630 +-0.270290 -0.178830 0.202360 0.493810 + +# -90 +1.456010 2.743180 2.589450 3.046230 3.451510 +3.319160 3.052900 3.873720 2.420650 0.949100 +0.008370 -1.382980 -2.138930 -2.087380 -1.268300 +-0.494370 0.267580 0.908250 0.537520 0.306260 +0.069540 0.097460 0.263060 0.603220 + +# -75 +1.396790 3.349090 2.180920 2.942960 3.814070 +3.675800 3.555310 3.887290 2.101260 -0.190940 +-0.732240 -1.382040 -0.673880 -0.817390 -0.826980 +-0.111800 0.053710 0.296400 0.692240 0.428960 +-0.036100 -0.033820 -0.194300 0.400210 + +# -60 +0.246650 1.229980 1.716960 3.168570 4.208190 +4.366860 4.251080 3.348110 0.997540 -1.287540 +-1.179900 -0.684300 -0.853660 -1.158760 -0.347550 +0.114810 0.242800 0.322420 0.370140 -0.374950 +-0.676940 -1.323430 -1.366650 -0.218770 + +# -45 +-1.196730 0.078060 2.347410 4.211350 5.376000 +5.364940 4.355200 2.436510 0.408470 -0.590840 +-0.435960 -0.501210 -0.822230 -0.607210 0.057910 +0.246580 -0.070570 0.379430 0.247770 -0.571680 +-1.282910 -1.715770 -1.839820 -1.987110 + +# -30 +-1.174720 1.067030 4.180460 6.741610 6.070770 +4.781470 2.758340 1.295810 0.571150 -0.196480 +0.251860 -0.732140 1.289360 1.497590 1.890550 +2.198490 0.169290 0.534000 0.331780 -1.276320 +-2.550070 -3.312150 -3.136670 -2.642260 + +# -15 +0.293590 5.588070 3.732620 3.217620 3.272450 +2.492320 1.563700 1.356760 0.831410 0.630170 +1.591970 0.821920 0.486070 0.715760 0.996020 +1.591580 -0.367400 0.181770 -0.613920 -2.267900 +-3.516460 -3.597700 -3.043340 -1.765020 + +# 0 +2.832310 0.787990 0.323280 0.479230 0.628600 +0.976330 1.238750 1.671950 1.645480 2.520340 +1.606970 0.776350 0.119780 0.070390 0.121170 +-1.569230 -1.213010 -1.846360 -2.744510 -3.792530 +-3.934880 -3.615930 -2.675750 -0.924170 + +# 15 +-0.778340 -1.912680 -2.052140 -1.846280 -1.047430 +0.183400 1.682950 2.223500 1.358370 2.448660 +1.436920 0.678570 -0.237060 -0.535320 -0.790380 +-2.182580 -3.251140 -4.195110 -4.269270 -3.908210 +-3.455620 -2.773970 1.755370 0.313410 + +# 30 +-2.963810 -3.483730 -3.517080 -2.724860 -1.405510 +0.336200 1.428450 1.394630 0.970370 2.462720 +1.522430 0.553620 -0.407380 -1.482950 -3.613920 +-4.159810 -4.945580 -4.784040 -3.764540 -2.959140 +-1.963850 -1.071260 -1.599580 -2.445320 + +# 45 +-4.029070 -3.932660 -3.558480 -2.513980 -1.037320 +0.362000 0.814380 0.754110 0.502370 1.903420 +0.770220 -0.416420 -3.286310 -3.875270 -4.907800 +-5.704430 -5.645660 -4.396040 -2.865450 -2.368170 +-2.860490 -3.416560 -3.666490 -3.859070 + +# 60 +-3.338270 -2.960220 -2.311700 -1.272890 -0.246470 +0.722610 0.668070 0.438130 2.395330 1.632470 +-2.041450 -3.218100 -3.915080 -4.852510 -5.696500 +-6.314370 -5.683690 -4.170620 -3.141000 -3.508820 +-3.756430 -3.640810 -3.640430 -3.550690 + +# 75 +-2.244860 -1.632100 -1.000640 -0.170440 0.526440 +0.823710 0.517140 -0.013120 -0.370910 -1.213720 +-2.305650 -3.420580 -4.484960 -5.693140 -6.199150 +-6.253870 -5.211310 -4.174380 -3.685150 -4.151360 +-4.161970 -3.725150 -3.715310 -2.606760 + +# 90 +-1.720840 -1.177830 -0.428430 0.277730 0.807900 +0.803260 0.482510 -0.336900 -0.786270 -1.774070 +-2.793220 -3.828560 -5.211800 -6.636850 -6.989940 +-6.108800 -5.452410 -3.911450 -4.321000 -4.587240 +-4.102610 -3.772820 -3.157300 -2.648390 + +# 105 +-1.850640 -1.092420 -0.445020 0.128490 1.005520 +0.884820 0.485850 -0.218470 -0.857670 -1.682330 +-3.014400 -4.481110 -6.053510 -6.865400 -6.871130 +-5.728240 -3.912230 -4.802110 -5.034640 -4.715990 +-4.601080 -4.086220 -3.274630 -2.410940 + +# 120 +-1.969230 -1.116650 -0.540250 -0.150330 0.763520 +1.038890 0.758480 0.313530 -0.333050 -1.872770 +-3.366270 -5.008260 -6.124810 -7.034830 -6.724320 +-3.700200 -4.510620 -5.185650 -5.361620 -4.847490 +-4.444320 -4.004260 -3.415720 -2.751230 + +# 135 +-2.111250 -1.168960 -0.322790 -0.006920 0.316660 +1.086270 0.939170 0.625340 -0.166360 -1.830310 +-3.469470 -4.946030 -6.112560 -1.915580 -4.047310 +-4.996740 -4.996730 -4.842690 -4.886620 -4.300540 +-4.494620 -4.442210 -4.163570 -3.183510 + +# 150 +-1.757590 -0.403620 0.023920 0.362390 0.634520 +1.264920 1.361360 0.948420 -0.073680 -1.483560 +-3.152820 1.835120 -1.762860 -5.093660 -5.744830 +-5.390070 -4.783930 -4.190630 -4.115420 -4.042280 +-4.125570 -4.028550 -4.026100 -2.937910 + +# 165 +-0.810590 -0.071500 0.378890 0.543310 1.277880 +1.641310 1.698840 1.519950 0.631950 -1.088670 +-2.736530 -0.735240 -4.563830 -6.408350 -5.889450 +-5.141750 -4.194970 -3.666490 -3.843450 -3.818830 +-3.826180 -3.596820 -2.994790 -2.231020 + +# alanine before proline map, type 2 + +# -180 +0.126790 0.768700 0.971260 1.250970 2.121010 +2.695430 2.064440 1.764790 0.755870 -0.713470 +0.976130 -2.475520 -5.455650 -5.096450 -5.305850 +-3.975630 -3.088580 -2.784200 -2.677120 -2.646060 +-2.335350 -2.010440 -1.608040 -0.482250 + +# -165 +-0.802290 1.377090 1.577020 1.872290 2.398990 +2.461630 2.333840 1.904070 1.061460 0.518400 +-0.116320 -3.575440 -5.284480 -5.160310 -4.196010 +-3.276210 -2.715340 -1.806200 -1.101780 -1.210320 +-1.008810 -0.637100 -1.603360 -1.776870 + +# -150 +-0.634810 1.156210 1.624350 2.047200 2.653910 +2.691410 2.296420 1.960450 1.324930 2.038290 +-1.151510 -3.148610 -4.058280 -4.531850 -3.796370 +-2.572090 -1.727250 -0.961410 -0.282910 -0.479120 +-1.039340 -1.618060 -1.725460 -1.376360 + +# -135 +0.214000 1.521370 1.977440 2.377950 2.929470 +2.893410 2.435810 2.162970 1.761500 1.190090 +-1.218610 -2.108900 -2.976100 -3.405340 -2.768440 +-1.836030 -0.957950 0.021790 -0.032760 -0.665880 +-1.321170 -1.212320 -0.893170 -0.897040 + +# -120 +0.873950 1.959160 2.508990 2.841100 3.698960 +3.309330 2.614300 2.481720 2.694660 1.082440 +-0.398320 -1.761800 -2.945110 -3.294690 -2.308300 +-0.855480 -0.087320 0.439040 0.691880 -0.586330 +-1.027210 -0.976640 -0.467580 0.104020 + +# -105 +1.767380 2.286650 2.818030 3.065500 3.370620 +3.397440 2.730310 2.878790 2.542010 1.545240 +-0.092150 -1.694440 -2.812310 -2.802430 -1.856360 +-0.306240 -0.122440 0.444680 0.810150 -0.058630 +-0.270290 -0.178830 0.202360 0.493810 + +# -90 +1.456010 2.743180 2.589450 3.046230 3.451510 +3.319160 3.052900 3.873720 2.420650 0.949100 +0.008370 -1.382980 -2.138930 -2.087380 -1.268300 +-0.494370 0.267580 0.908250 0.537520 0.306260 +0.069540 0.097460 0.263060 0.603220 + +# -75 +1.396790 3.349090 2.180920 2.942960 3.814070 +3.675800 3.555310 3.887290 2.101260 -0.190940 +-0.732240 -1.382040 -0.673880 -0.817390 -0.826980 +-0.111800 0.053710 0.296400 0.692240 0.428960 +-0.036100 -0.033820 -0.194300 0.400210 + +# -60 +0.246650 1.229980 1.716960 3.168570 4.208190 +4.366860 4.251080 3.348110 0.997540 -1.287540 +-1.179900 -0.684300 -0.853660 -1.158760 -0.347550 +0.114810 0.242800 0.322420 0.370140 -0.374950 +-0.676940 -1.323430 -1.366650 -0.218770 + +# -45 +-1.196730 0.078060 2.347410 4.211350 5.376000 +5.364940 4.355200 2.436510 0.408470 -0.590840 +-0.435960 -0.501210 -0.822230 -0.607210 0.057910 +0.246580 -0.070570 0.379430 0.247770 -0.571680 +-1.282910 -1.715770 -1.839820 -1.987110 + +# -30 +-1.174720 1.067030 4.180460 6.741610 6.070770 +4.781470 2.758340 1.295810 0.571150 -0.196480 +0.251860 -0.732140 1.289360 1.497590 1.890550 +2.198490 0.169290 0.534000 0.331780 -1.276320 +-2.550070 -3.312150 -3.136670 -2.642260 + +# -15 +0.293590 5.588070 3.732620 3.217620 3.272450 +2.492320 1.563700 1.356760 0.831410 0.630170 +1.591970 0.821920 0.486070 0.715760 0.996020 +1.591580 -0.367400 0.181770 -0.613920 -2.267900 +-3.516460 -3.597700 -3.043340 -1.765020 + +# 0 +2.832310 0.787990 0.323280 0.479230 0.628600 +0.976330 1.238750 1.671950 1.645480 2.520340 +1.606970 0.776350 0.119780 0.070390 0.121170 +-1.569230 -1.213010 -1.846360 -2.744510 -3.792530 +-3.934880 -3.615930 -2.675750 -0.924170 + +# 15 +-0.778340 -1.912680 -2.052140 -1.846280 -1.047430 +0.183400 1.682950 2.223500 1.358370 2.448660 +1.436920 0.678570 -0.237060 -0.535320 -0.790380 +-2.182580 -3.251140 -4.195110 -4.269270 -3.908210 +-3.455620 -2.773970 1.755370 0.313410 + +# 30 +-2.963810 -3.483730 -3.517080 -2.724860 -1.405510 +0.336200 1.428450 1.394630 0.970370 2.462720 +1.522430 0.553620 -0.407380 -1.482950 -3.613920 +-4.159810 -4.945580 -4.784040 -3.764540 -2.959140 +-1.963850 -1.071260 -1.599580 -2.445320 + +# 45 +-4.029070 -3.932660 -3.558480 -2.513980 -1.037320 +0.362000 0.814380 0.754110 0.502370 1.903420 +0.770220 -0.416420 -3.286310 -3.875270 -4.907800 +-5.704430 -5.645660 -4.396040 -2.865450 -2.368170 +-2.860490 -3.416560 -3.666490 -3.859070 + +# 60 +-3.338270 -2.960220 -2.311700 -1.272890 -0.246470 +0.722610 0.668070 0.438130 2.395330 1.632470 +-2.041450 -3.218100 -3.915080 -4.852510 -5.696500 +-6.314370 -5.683690 -4.170620 -3.141000 -3.508820 +-3.756430 -3.640810 -3.640430 -3.550690 + +# 75 +-2.244860 -1.632100 -1.000640 -0.170440 0.526440 +0.823710 0.517140 -0.013120 -0.370910 -1.213720 +-2.305650 -3.420580 -4.484960 -5.693140 -6.199150 +-6.253870 -5.211310 -4.174380 -3.685150 -4.151360 +-4.161970 -3.725150 -3.715310 -2.606760 + +# 90 +-1.720840 -1.177830 -0.428430 0.277730 0.807900 +0.803260 0.482510 -0.336900 -0.786270 -1.774070 +-2.793220 -3.828560 -5.211800 -6.636850 -6.989940 +-6.108800 -5.452410 -3.911450 -4.321000 -4.587240 +-4.102610 -3.772820 -3.157300 -2.648390 + +# 105 +-1.850640 -1.092420 -0.445020 0.128490 1.005520 +0.884820 0.485850 -0.218470 -0.857670 -1.682330 +-3.014400 -4.481110 -6.053510 -6.865400 -6.871130 +-5.728240 -3.912230 -4.802110 -5.034640 -4.715990 +-4.601080 -4.086220 -3.274630 -2.410940 + +# 120 +-1.969230 -1.116650 -0.540250 -0.150330 0.763520 +1.038890 0.758480 0.313530 -0.333050 -1.872770 +-3.366270 -5.008260 -6.124810 -7.034830 -6.724320 +-3.700200 -4.510620 -5.185650 -5.361620 -4.847490 +-4.444320 -4.004260 -3.415720 -2.751230 + +# 135 +-2.111250 -1.168960 -0.322790 -0.006920 0.316660 +1.086270 0.939170 0.625340 -0.166360 -1.830310 +-3.469470 -4.946030 -6.112560 -1.915580 -4.047310 +-4.996740 -4.996730 -4.842690 -4.886620 -4.300540 +-4.494620 -4.442210 -4.163570 -3.183510 + +# 150 +-1.757590 -0.403620 0.023920 0.362390 0.634520 +1.264920 1.361360 0.948420 -0.073680 -1.483560 +-3.152820 1.835120 -1.762860 -5.093660 -5.744830 +-5.390070 -4.783930 -4.190630 -4.115420 -4.042280 +-4.125570 -4.028550 -4.026100 -2.937910 + +# 165 +-0.810590 -0.071500 0.378890 0.543310 1.277880 +1.641310 1.698840 1.519950 0.631950 -1.088670 +-2.736530 -0.735240 -4.563830 -6.408350 -5.889450 +-5.141750 -4.194970 -3.666490 -3.843450 -3.818830 +-3.826180 -3.596820 -2.994790 -2.231020 + +# proline, type 3 + +# -180.00 +-4.60660 -4.28920 -4.51420 -5.49210 -6.65460 +-7.18530 -7.63320 -8.41920 -9.12510 -8.97830 +-9.02750 -8.88890 -8.61060 -8.10240 -7.96680 +-7.98860 -7.96190 -7.56330 -7.73950 -8.31580 +-9.33380 -9.61880 -7.91860 -6.00570 + +# -165.00 +-3.21030 -2.90350 -3.04340 -3.83720 -4.92360 +-5.41790 -5.78320 -6.62490 -7.17820 -7.59500 +-7.60980 -7.39440 -6.68190 -6.24950 -5.83810 +-5.90510 -5.57900 -4.89950 -5.18840 -6.13180 +-6.93510 -7.49160 -5.84410 -4.48890 + +# -150.00 +-2.23750 -2.26520 -2.55120 -3.67610 -4.61620 +-5.01560 -5.36920 -5.93110 -6.16820 -6.48830 +-6.54220 -6.14310 -5.33350 -4.85090 -4.50790 +-4.44050 -4.34380 -4.13820 -4.35620 -5.24880 +-6.05280 -6.18610 -4.80390 -3.41730 + +# -135.00 +-2.12830 -2.02010 -2.90220 -3.88240 -4.63610 +-5.02000 -5.46090 -5.79830 -5.63220 -6.41490 +-6.20330 -5.57550 -4.58160 -4.15850 -4.00710 +-3.95020 -3.76880 -4.09690 -4.52680 -5.33410 +-5.97380 -5.68710 -4.70040 -3.33790 + +# -120.00 +-2.20610 -2.22370 -2.85110 -3.55370 -4.33320 +-4.57860 -4.95030 -5.23950 -5.37730 -6.14740 +-5.98660 -5.17760 -4.16090 -3.78450 -3.83210 +-3.86590 -4.21340 -4.26250 -4.16200 -4.62160 +-5.47440 -5.38280 -4.39100 -3.16760 + +# -105.00 +-1.32980 -1.44660 -2.00420 -2.84440 -3.31580 +-3.15410 -2.83060 -3.04980 -4.04810 -5.02030 +-4.84090 -3.88520 -2.79540 -2.49190 -2.65990 +-3.34670 -3.60580 -3.76970 -4.07430 -4.21930 +-4.66740 -4.23940 -3.37160 -2.10050 + +# -90.00 +-0.27400 -0.41780 -0.92340 -1.48470 -1.84370 +-1.42480 -1.14210 -1.83280 -2.97900 -3.57000 +-3.33310 -2.24200 -1.31350 -0.96640 0.38300 +-0.57540 -0.74080 -0.59130 -0.41960 -2.64120 +-3.25380 -2.85540 -2.09970 -0.95660 + +# -75.00 + 0.26880 0.18190 0.00300 -0.45390 -0.24210 + 0.21240 0.33320 -0.95520 -2.17670 -2.66490 +-2.24910 -1.35480 -0.35030 1.84000 0.81880 + 0.41530 0.58460 0.65350 0.54960 0.41300 + 0.23350 -2.13490 -1.54900 -0.38340 + +# -60.00 +-0.48820 -0.40730 -0.25210 -0.20490 0.53370 + 1.05990 0.72540 -0.68300 -2.10190 -2.83530 +-2.55500 -1.33940 -0.24710 1.63190 1.10600 + 1.16350 1.30720 1.49310 1.16160 0.99320 +-2.33360 -2.59380 -2.01110 -1.46840 + +# -45.00 +-2.41010 -1.30740 -0.85960 -0.06830 0.20500 + 0.00650 -0.59750 -2.39820 -3.98270 -4.46610 +-3.51230 -1.99330 0.96050 0.80870 0.64450 + 0.82590 1.00550 1.28110 0.96670 -2.71140 +-4.18330 -4.83370 -4.49950 -3.47920 + +# -30.00 +-3.46300 -2.06680 -0.84260 -0.70440 -1.85140 +-2.59990 -3.62600 -5.08050 -6.35480 -6.08340 +-4.78650 -1.24210 -0.81660 -0.60170 -0.50040 +-0.41020 -0.28450 -0.14620 -3.14950 -5.16280 +-6.81940 -7.52060 -7.12480 -5.49420 + +# -15.00 +-2.70550 -1.11410 -0.82120 -3.80470 -4.50780 +-4.64960 -5.78530 -7.04450 -7.30810 -6.56450 +-3.47860 -2.68630 -1.78120 -1.48400 -1.45930 +-1.43060 -5.75320 -2.72500 -6.45790 -8.16530 +-9.86330 -9.80710 -8.62670 -6.49690 + +# 0.00 +-2.73190 -3.83060 -5.26960 -5.91610 -6.71690 +-7.05420 -7.77520 -7.81220 -7.54160 -4.98400 +-5.05070 -4.28910 -3.78740 -3.64600 -3.91990 +-3.55740 -3.03470 -5.78910 -7.02720 -8.04090 +-8.60950 -8.27740 -7.10120 -5.12040 + +# 15.00 +-3.17420 -4.66830 -6.17790 -7.26760 -8.22560 +-8.32120 -7.64880 -7.00240 -4.77210 -5.65380 +-5.36790 -4.43640 -4.37720 -4.45280 -4.47980 +-4.07260 -6.75420 -7.56910 -7.95830 -7.87760 +-7.88160 -7.29010 -7.15600 -5.19260 + +# 30.00 +-2.96240 -4.87720 -6.64930 -7.95180 -8.47610 +-7.49950 -6.26260 -6.35040 -4.34190 -4.49790 +-4.27480 -1.75620 -3.89990 -3.81400 -3.99610 +-6.66430 -7.66400 -7.80660 -7.21980 -6.83750 +-6.88220 -7.10200 -5.10200 -2.38920 + +# 45 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 60 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 75 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 90 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 105 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 120 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 135 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 150 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 165 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 2 adjacent prolines, type 4 + +# -180.00 +-4.60660 -4.28920 -4.51420 -5.49210 -6.65460 +-7.18530 -7.63320 -8.41920 -9.12510 -8.97830 +-9.02750 -8.88890 -8.61060 -8.10240 -7.96680 +-7.98860 -7.96190 -7.56330 -7.73950 -8.31580 +-9.33380 -9.61880 -7.91860 -6.00570 + +# -165.00 +-3.21030 -2.90350 -3.04340 -3.83720 -4.92360 +-5.41790 -5.78320 -6.62490 -7.17820 -7.59500 +-7.60980 -7.39440 -6.68190 -6.24950 -5.83810 +-5.90510 -5.57900 -4.89950 -5.18840 -6.13180 +-6.93510 -7.49160 -5.84410 -4.48890 + +# -150.00 +-2.23750 -2.26520 -2.55120 -3.67610 -4.61620 +-5.01560 -5.36920 -5.93110 -6.16820 -6.48830 +-6.54220 -6.14310 -5.33350 -4.85090 -4.50790 +-4.44050 -4.34380 -4.13820 -4.35620 -5.24880 +-6.05280 -6.18610 -4.80390 -3.41730 + +# -135.00 +-2.12830 -2.02010 -2.90220 -3.88240 -4.63610 +-5.02000 -5.46090 -5.79830 -5.63220 -6.41490 +-6.20330 -5.57550 -4.58160 -4.15850 -4.00710 +-3.95020 -3.76880 -4.09690 -4.52680 -5.33410 +-5.97380 -5.68710 -4.70040 -3.33790 + +# -120.00 +-2.20610 -2.22370 -2.85110 -3.55370 -4.33320 +-4.57860 -4.95030 -5.23950 -5.37730 -6.14740 +-5.98660 -5.17760 -4.16090 -3.78450 -3.83210 +-3.86590 -4.21340 -4.26250 -4.16200 -4.62160 +-5.47440 -5.38280 -4.39100 -3.16760 + +# -105.00 +-1.32980 -1.44660 -2.00420 -2.84440 -3.31580 +-3.15410 -2.83060 -3.04980 -4.04810 -5.02030 +-4.84090 -3.88520 -2.79540 -2.49190 -2.65990 +-3.34670 -3.60580 -3.76970 -4.07430 -4.21930 +-4.66740 -4.23940 -3.37160 -2.10050 + +# -90.00 +-0.27400 -0.41780 -0.92340 -1.48470 -1.84370 +-1.42480 -1.14210 -1.83280 -2.97900 -3.57000 +-3.33310 -2.24200 -1.31350 -0.96640 0.38300 +-0.57540 -0.74080 -0.59130 -0.41960 -2.64120 +-3.25380 -2.85540 -2.09970 -0.95660 + +# -75.00 + 0.26880 0.18190 0.00300 -0.45390 -0.24210 + 0.21240 0.33320 -0.95520 -2.17670 -2.66490 +-2.24910 -1.35480 -0.35030 1.84000 0.81880 + 0.41530 0.58460 0.65350 0.54960 0.41300 + 0.23350 -2.13490 -1.54900 -0.38340 + +# -60.00 +-0.48820 -0.40730 -0.25210 -0.20490 0.53370 + 1.05990 0.72540 -0.68300 -2.10190 -2.83530 +-2.55500 -1.33940 -0.24710 1.63190 1.10600 + 1.16350 1.30720 1.49310 1.16160 0.99320 +-2.33360 -2.59380 -2.01110 -1.46840 + +# -45.00 +-2.41010 -1.30740 -0.85960 -0.06830 0.20500 + 0.00650 -0.59750 -2.39820 -3.98270 -4.46610 +-3.51230 -1.99330 0.96050 0.80870 0.64450 + 0.82590 1.00550 1.28110 0.96670 -2.71140 +-4.18330 -4.83370 -4.49950 -3.47920 + +# -30.00 +-3.46300 -2.06680 -0.84260 -0.70440 -1.85140 +-2.59990 -3.62600 -5.08050 -6.35480 -6.08340 +-4.78650 -1.24210 -0.81660 -0.60170 -0.50040 +-0.41020 -0.28450 -0.14620 -3.14950 -5.16280 +-6.81940 -7.52060 -7.12480 -5.49420 + +# -15.00 +-2.70550 -1.11410 -0.82120 -3.80470 -4.50780 +-4.64960 -5.78530 -7.04450 -7.30810 -6.56450 +-3.47860 -2.68630 -1.78120 -1.48400 -1.45930 +-1.43060 -5.75320 -2.72500 -6.45790 -8.16530 +-9.86330 -9.80710 -8.62670 -6.49690 + +# 0.00 +-2.73190 -3.83060 -5.26960 -5.91610 -6.71690 +-7.05420 -7.77520 -7.81220 -7.54160 -4.98400 +-5.05070 -4.28910 -3.78740 -3.64600 -3.91990 +-3.55740 -3.03470 -5.78910 -7.02720 -8.04090 +-8.60950 -8.27740 -7.10120 -5.12040 + +# 15.00 +-3.17420 -4.66830 -6.17790 -7.26760 -8.22560 +-8.32120 -7.64880 -7.00240 -4.77210 -5.65380 +-5.36790 -4.43640 -4.37720 -4.45280 -4.47980 +-4.07260 -6.75420 -7.56910 -7.95830 -7.87760 +-7.88160 -7.29010 -7.15600 -5.19260 + +# 30.00 +-2.96240 -4.87720 -6.64930 -7.95180 -8.47610 +-7.49950 -6.26260 -6.35040 -4.34190 -4.49790 +-4.27480 -1.75620 -3.89990 -3.81400 -3.99610 +-6.66430 -7.66400 -7.80660 -7.21980 -6.83750 +-6.88220 -7.10200 -5.10200 -2.38920 + +# 45 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 60 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 75 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 90 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 105 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 120 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 135 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 150 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# 165 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 0.000000 +0.000000 0.000000 0.000000 0.000000 + +# glycine map, type 5 + +# -180 +-0.549160 -0.535500 -0.588110 -0.754620 -0.679290 +-0.038150 0.298460 0.326040 -0.375610 -1.704360 +-3.061280 -3.956460 -3.576280 -1.038930 2.012450 +-1.714610 -0.377660 0.317310 0.294580 -0.042920 +-0.676620 -0.744600 -0.586590 -0.554770 + +# -165 +-0.709450 -0.896700 -0.990850 -1.319240 -0.940260 +-0.126160 0.329180 0.258100 -0.534910 -1.715700 +-2.780320 -3.153350 -1.636020 1.822690 -2.675640 +-1.810120 -0.410680 0.180860 0.196710 -0.000430 +-0.271890 -0.462500 -0.348750 -0.477660 + +# -150 +-1.224850 -1.482430 -1.665900 -1.656770 -1.119780 +-1.642540 -0.054220 -0.290670 -0.887080 -1.626260 +-2.165440 -1.546500 0.753400 -2.949180 -2.225630 +-1.664160 -0.628990 0.000490 0.033160 -0.092820 +-0.339050 -0.563330 -0.794980 -0.710760 + +# -135 +-1.787640 -2.117750 -2.143020 -1.803720 -1.567160 +-0.886880 -0.801350 -0.851590 -1.020630 -1.337360 +-1.062570 0.338010 -4.372310 -2.435890 -2.220710 +-1.718060 -0.758950 -0.207560 0.100910 -0.055650 +-0.288370 -0.880610 -1.267450 -1.465530 + +# -120 +-2.348270 -2.593790 -2.596140 -2.364070 -1.970070 +-1.705860 -1.435540 -1.289220 -1.358170 -0.975570 +-3.514390 -4.283210 -3.975820 -3.215190 -2.394430 +-1.455320 -0.553910 -0.158900 -0.173830 -0.297950 +-0.661220 -1.068330 -1.601800 -1.914850 + +# -105 +-2.788800 -3.079570 -3.178150 -3.013710 -2.626630 +-2.266680 -1.951490 -1.681850 -1.195390 -2.567680 +-3.632800 -4.748210 -4.662850 -4.255190 -2.776760 +-1.695490 -0.893140 -0.633810 -0.467320 -0.540540 +-0.950190 -1.401500 -1.959970 -2.412680 + +# -90 +-3.857170 -3.713610 -3.902110 -3.611370 -3.040850 +-2.406460 -1.975250 -1.452040 -0.971860 -2.808170 +-4.181160 -4.981430 -5.446890 -4.359900 -2.864390 +-1.898510 -1.139090 -0.971340 -1.065550 -1.020680 +-1.141350 -1.794480 -2.420970 -2.939990 + +# -75 +-4.987770 -4.995210 -4.485310 -3.892550 -3.228630 +-2.345360 -1.664160 -1.105500 -1.945510 -3.715530 +-4.492140 -5.536170 -5.708500 -3.675410 -2.986660 +-1.859410 -0.756620 -1.269930 -1.312730 -1.607440 +-1.892510 -2.659400 -3.347950 -3.970600 + +# -60 +-6.183650 -5.456080 -4.878940 -4.000820 -2.683230 +-2.067520 -1.094850 -1.119790 -2.962910 -3.687830 +-4.993340 -4.666260 -3.796280 -3.374140 -2.495430 +-1.453990 -0.877560 -1.002930 -1.337310 -2.431360 +-2.948140 -4.008100 -4.821040 -5.565810 + +# -45 +-6.755760 -5.850030 -4.362190 -2.714090 -1.708710 +-0.526660 -0.536700 -2.037170 -3.892650 -4.558570 +-4.237410 -3.735160 -3.688580 -3.009910 -2.112940 +-1.455400 -0.925490 -1.121840 -1.561900 -2.751370 +-4.094860 -5.207530 -6.128530 -6.613030 + +# -30 +-5.716250 -4.434060 -2.788600 -0.974400 -0.729200 +-0.904940 -1.833540 -3.017700 -3.313450 -3.336010 +-3.181640 -3.594720 -1.231370 -0.603790 0.128810 +-1.222610 -0.909150 -0.837700 -1.346820 -3.040880 +-4.731110 -5.844860 -6.428460 -6.424880 + +# -15 +-3.991110 -2.046000 0.082550 -2.676110 -2.828500 +-2.596640 -2.843330 -3.011480 -2.312640 -2.405980 +-3.086210 -1.164620 -1.231660 -0.871900 -0.348980 +-1.735900 -0.914150 -0.484520 -1.818040 -3.602550 +-5.330320 -5.992270 -5.588080 -5.408360 + +# 0 +-1.147060 -3.317730 -4.305100 -4.615200 -4.533780 +-3.622950 -2.832800 -1.872810 -1.144300 -1.994070 +-0.741980 -1.115010 -1.229250 -1.103680 -0.742430 +-1.973970 -1.070020 -1.802220 -2.712770 -3.624130 +-4.537100 -4.619970 -4.310890 -3.318290 + +# 15 +-3.997710 -5.408360 -5.588080 -5.992270 -5.330320 +-3.602550 -1.818040 -0.484520 -0.914150 -1.735900 +-0.348980 -0.871900 -1.231660 -1.164620 -3.086210 +-2.405980 -2.312640 -3.011480 -2.843330 -2.596640 +-2.828500 -2.676110 0.082550 -2.046000 + +# 30 +-5.710850 -6.424880 -6.428460 -5.844860 -4.731110 +-3.040880 -1.346820 -0.837700 -0.909150 -1.222610 +0.128810 -0.603790 -1.231370 -3.594720 -3.181640 +-3.336010 -3.313450 -3.017700 -1.833540 -0.904940 +-0.729200 -0.974400 -2.788600 -4.434060 + +# 45 +-6.754940 -6.613030 -6.128530 -5.207530 -4.094860 +-2.751370 -1.561900 -1.121840 -0.925490 -1.455400 +-2.112940 -3.009910 -3.688580 -3.735160 -4.237410 +-4.558570 -3.892650 -2.037170 -0.536700 -0.526660 +-1.708710 -2.714090 -4.362190 -5.850030 + +# 60 +-6.188070 -5.565810 -4.821040 -4.008100 -2.948140 +-2.431360 -1.337310 -1.002930 -0.877560 -1.453990 +-2.495430 -3.374140 -3.796280 -4.666260 -4.993340 +-3.687830 -2.962910 -1.119790 -1.094850 -2.067520 +-2.683230 -4.000820 -4.878940 -5.456080 + +# 75 +-4.986080 -3.970600 -3.347950 -2.659400 -1.892510 +-1.607440 -1.312730 -1.269930 -0.756620 -1.859410 +-2.986660 -3.675410 -5.708500 -5.536170 -4.492140 +-3.715530 -1.945510 -1.105500 -1.664160 -2.345360 +-3.228630 -3.892550 -4.485310 -4.995210 + +# 90 +-3.879190 -2.939990 -2.420970 -1.794480 -1.141350 +-1.020680 -1.065550 -0.971340 -1.139090 -1.898510 +-2.864390 -4.359900 -5.446890 -4.981430 -4.181160 +-2.808170 -0.971860 -1.452040 -1.975250 -2.406460 +-3.040850 -3.611370 -3.902110 -3.713610 + +# 105 +-2.793280 -2.412680 -1.959970 -1.401500 -0.950190 +-0.540540 -0.467320 -0.633810 -0.893140 -1.695490 +-2.776760 -4.255190 -4.662850 -4.448210 -3.332800 +-2.567680 -1.195390 -1.681850 -1.951490 -2.266680 +-2.626630 -3.013710 -3.178150 -3.079570 + +# 120 +-2.330190 -1.914850 -1.601800 -1.068330 -0.661220 +-0.297950 -0.173830 -0.158900 -0.553910 -1.455320 +-2.394430 -3.215190 -3.975820 -3.783210 -3.014390 +-0.975570 -1.358170 -1.289220 -1.435540 -1.705860 +-1.970070 -2.364070 -2.596140 -2.593790 + +# 135 +-1.796120 -1.465530 -1.267450 -0.880610 -0.288370 +-0.055650 0.100910 -0.207560 -0.758950 -1.718060 +-2.220710 -2.435890 -4.372310 0.338010 -1.062570 +-1.337360 -1.020630 -0.851590 -0.801350 -0.886880 +-1.567160 -1.803720 -2.143020 -2.117750 + +# 150 +-1.263610 -0.710760 -0.794980 -0.563330 -0.339050 +-0.092820 0.033160 0.000490 -0.628990 -1.664160 +-2.225630 -2.949180 0.753400 -1.546500 -2.165440 +-1.626260 -0.887080 -0.290670 -0.054220 -1.642540 +-1.119780 -1.656770 -1.665900 -1.482430 + +# 165 +-0.684660 -0.477660 -0.348750 -0.462500 -0.271890 +-0.000430 0.196710 0.180860 -0.410680 -1.810120 +-2.675640 1.822690 -1.636020 -3.153350 -2.780320 +-1.715700 -0.534910 0.258100 0.329180 -0.126160 +-0.940260 -1.319240 -0.990850 -0.896700 + +# glycine before proline map, type 6 + +# -180 +-0.549160 -0.535500 -0.588110 -0.754620 -0.679290 +-0.038150 0.298460 0.326040 -0.375610 -1.704360 +-3.061280 -3.956460 -3.576280 -1.038930 2.012450 +-1.714610 -0.377660 0.317310 0.294580 -0.042920 +-0.676620 -0.744600 -0.586590 -0.554770 + +# -165 +-0.709450 -0.896700 -0.990850 -1.319240 -0.940260 +-0.126160 0.329180 0.258100 -0.534910 -1.715700 +-2.780320 -3.153350 -1.636020 1.822690 -2.675640 +-1.810120 -0.410680 0.180860 0.196710 -0.000430 +-0.271890 -0.462500 -0.348750 -0.477660 + +# -150 +-1.224850 -1.482430 -1.665900 -1.656770 -1.119780 +-1.642540 -0.054220 -0.290670 -0.887080 -1.626260 +-2.165440 -1.546500 0.753400 -2.949180 -2.225630 +-1.664160 -0.628990 0.000490 0.033160 -0.092820 +-0.339050 -0.563330 -0.794980 -0.710760 + +# -135 +-1.787640 -2.117750 -2.143020 -1.803720 -1.567160 +-0.886880 -0.801350 -0.851590 -1.020630 -1.337360 +-1.062570 0.338010 -4.372310 -2.435890 -2.220710 +-1.718060 -0.758950 -0.207560 0.100910 -0.055650 +-0.288370 -0.880610 -1.267450 -1.465530 + +# -120 +-2.348270 -2.593790 -2.596140 -2.364070 -1.970070 +-1.705860 -1.435540 -1.289220 -1.358170 -0.975570 +-3.514390 -4.283210 -3.975820 -3.215190 -2.394430 +-1.455320 -0.553910 -0.158900 -0.173830 -0.297950 +-0.661220 -1.068330 -1.601800 -1.914850 + +# -105 +-2.788800 -3.079570 -3.178150 -3.013710 -2.626630 +-2.266680 -1.951490 -1.681850 -1.195390 -2.567680 +-3.632800 -4.748210 -4.662850 -4.255190 -2.776760 +-1.695490 -0.893140 -0.633810 -0.467320 -0.540540 +-0.950190 -1.401500 -1.959970 -2.412680 + +# -90 +-3.857170 -3.713610 -3.902110 -3.611370 -3.040850 +-2.406460 -1.975250 -1.452040 -0.971860 -2.808170 +-4.181160 -4.981430 -5.446890 -4.359900 -2.864390 +-1.898510 -1.139090 -0.971340 -1.065550 -1.020680 +-1.141350 -1.794480 -2.420970 -2.939990 + +# -75 +-4.987770 -4.995210 -4.485310 -3.892550 -3.228630 +-2.345360 -1.664160 -1.105500 -1.945510 -3.715530 +-4.492140 -5.536170 -5.708500 -3.675410 -2.986660 +-1.859410 -0.756620 -1.269930 -1.312730 -1.607440 +-1.892510 -2.659400 -3.347950 -3.970600 + +# -60 +-6.183650 -5.456080 -4.878940 -4.000820 -2.683230 +-2.067520 -1.094850 -1.119790 -2.962910 -3.687830 +-4.993340 -4.666260 -3.796280 -3.374140 -2.495430 +-1.453990 -0.877560 -1.002930 -1.337310 -2.431360 +-2.948140 -4.008100 -4.821040 -5.565810 + +# -45 +-6.755760 -5.850030 -4.362190 -2.714090 -1.708710 +-0.526660 -0.536700 -2.037170 -3.892650 -4.558570 +-4.237410 -3.735160 -3.688580 -3.009910 -2.112940 +-1.455400 -0.925490 -1.121840 -1.561900 -2.751370 +-4.094860 -5.207530 -6.128530 -6.613030 + +# -30 +-5.716250 -4.434060 -2.788600 -0.974400 -0.729200 +-0.904940 -1.833540 -3.017700 -3.313450 -3.336010 +-3.181640 -3.594720 -1.231370 -0.603790 0.128810 +-1.222610 -0.909150 -0.837700 -1.346820 -3.040880 +-4.731110 -5.844860 -6.428460 -6.424880 + +# -15 +-3.991110 -2.046000 0.082550 -2.676110 -2.828500 +-2.596640 -2.843330 -3.011480 -2.312640 -2.405980 +-3.086210 -1.164620 -1.231660 -0.871900 -0.348980 +-1.735900 -0.914150 -0.484520 -1.818040 -3.602550 +-5.330320 -5.992270 -5.588080 -5.408360 + +# 0 +-1.147060 -3.317730 -4.305100 -4.615200 -4.533780 +-3.622950 -2.832800 -1.872810 -1.144300 -1.994070 +-0.741980 -1.115010 -1.229250 -1.103680 -0.742430 +-1.973970 -1.070020 -1.802220 -2.712770 -3.624130 +-4.537100 -4.619970 -4.310890 -3.318290 + +# 15 +-3.997710 -5.408360 -5.588080 -5.992270 -5.330320 +-3.602550 -1.818040 -0.484520 -0.914150 -1.735900 +-0.348980 -0.871900 -1.231660 -1.164620 -3.086210 +-2.405980 -2.312640 -3.011480 -2.843330 -2.596640 +-2.828500 -2.676110 0.082550 -2.046000 + +# 30 +-5.710850 -6.424880 -6.428460 -5.844860 -4.731110 +-3.040880 -1.346820 -0.837700 -0.909150 -1.222610 +0.128810 -0.603790 -1.231370 -3.594720 -3.181640 +-3.336010 -3.313450 -3.017700 -1.833540 -0.904940 +-0.729200 -0.974400 -2.788600 -4.434060 + +# 45 +-6.754940 -6.613030 -6.128530 -5.207530 -4.094860 +-2.751370 -1.561900 -1.121840 -0.925490 -1.455400 +-2.112940 -3.009910 -3.688580 -3.735160 -4.237410 +-4.558570 -3.892650 -2.037170 -0.536700 -0.526660 +-1.708710 -2.714090 -4.362190 -5.850030 + +# 60 +-6.188070 -5.565810 -4.821040 -4.008100 -2.948140 +-2.431360 -1.337310 -1.002930 -0.877560 -1.453990 +-2.495430 -3.374140 -3.796280 -4.666260 -4.993340 +-3.687830 -2.962910 -1.119790 -1.094850 -2.067520 +-2.683230 -4.000820 -4.878940 -5.456080 + +# 75 +-4.986080 -3.970600 -3.347950 -2.659400 -1.892510 +-1.607440 -1.312730 -1.269930 -0.756620 -1.859410 +-2.986660 -3.675410 -5.708500 -5.536170 -4.492140 +-3.715530 -1.945510 -1.105500 -1.664160 -2.345360 +-3.228630 -3.892550 -4.485310 -4.995210 + +# 90 +-3.879190 -2.939990 -2.420970 -1.794480 -1.141350 +-1.020680 -1.065550 -0.971340 -1.139090 -1.898510 +-2.864390 -4.359900 -5.446890 -4.981430 -4.181160 +-2.808170 -0.971860 -1.452040 -1.975250 -2.406460 +-3.040850 -3.611370 -3.902110 -3.713610 + +# 105 +-2.793280 -2.412680 -1.959970 -1.401500 -0.950190 +-0.540540 -0.467320 -0.633810 -0.893140 -1.695490 +-2.776760 -4.255190 -4.662850 -4.448210 -3.332800 +-2.567680 -1.195390 -1.681850 -1.951490 -2.266680 +-2.626630 -3.013710 -3.178150 -3.079570 + +# 120 +-2.330190 -1.914850 -1.601800 -1.068330 -0.661220 +-0.297950 -0.173830 -0.158900 -0.553910 -1.455320 +-2.394430 -3.215190 -3.975820 -3.783210 -3.014390 +-0.975570 -1.358170 -1.289220 -1.435540 -1.705860 +-1.970070 -2.364070 -2.596140 -2.593790 + +# 135 +-1.796120 -1.465530 -1.267450 -0.880610 -0.288370 +-0.055650 0.100910 -0.207560 -0.758950 -1.718060 +-2.220710 -2.435890 -4.372310 0.338010 -1.062570 +-1.337360 -1.020630 -0.851590 -0.801350 -0.886880 +-1.567160 -1.803720 -2.143020 -2.117750 + +# 150 +-1.263610 -0.710760 -0.794980 -0.563330 -0.339050 +-0.092820 0.033160 0.000490 -0.628990 -1.664160 +-2.225630 -2.949180 0.753400 -1.546500 -2.165440 +-1.626260 -0.887080 -0.290670 -0.054220 -1.642540 +-1.119780 -1.656770 -1.665900 -1.482430 + +# 165 +-0.684660 -0.477660 -0.348750 -0.462500 -0.271890 +-0.000430 0.196710 0.180860 -0.410680 -1.810120 +-2.675640 1.822690 -1.636020 -3.153350 -2.780320 +-1.715700 -0.534910 0.258100 0.329180 -0.126160 +-0.940260 -1.319240 -0.990850 -0.896700 + diff --git a/examples/cmap/gagg.data b/examples/cmap/gagg.data new file mode 100644 index 0000000000..fd2aa193c4 --- /dev/null +++ b/examples/cmap/gagg.data @@ -0,0 +1,380 @@ +Created by charmm2lammps v1.8.2.6 beta on Sun Mar 20 00:26:35 EDT 2016 + + 34 atoms + 33 bonds + 57 angles + 75 dihedrals + 7 impropers + 2 crossterms + + 13 atom types + 15 bond types + 30 angle types + 42 dihedral types + 5 improper types + + -34.414709 45.585291 xlo xhi + -36.134827 43.865173 ylo yhi + -39.349142 40.650858 zlo zhi + +Masses + + 1 1.008 # H + 2 1.008 # HC + 3 1.008 # HA + 4 1.008 # HB + 5 12.011 # C + 6 12.011 # CT1 + 7 12.011 # CT2 + 8 12.011 # CT3 + 9 12.011 # CC + 10 14.007 # NH1 + 11 14.007 # NH3 + 12 15.999 # O + 13 15.999 # OC + +Pair Coeffs + + 1 0.046 0.400013524445 0.046 0.400013524445 # H + 2 0.046 0.400013524445 0.046 0.400013524445 # HC + 3 0.022 2.35197261589 0.022 2.35197261589 # HA + 4 0.022 2.35197261589 0.022 2.35197261589 # HB + 5 0.11 3.56359487256 0.11 3.56359487256 # C + 6 0.02 4.05358916754 0.01 3.38541512893 # CT1 + 7 0.055 3.87540942391 0.01 3.38541512893 # CT2 + 8 0.08 3.67050271874 0.01 3.38541512893 # CT3 + 9 0.07 3.56359487256 0.07 3.56359487256 # CC + 10 0.2 3.29632525712 0.2 2.76178602624 # NH1 + 11 0.2 3.29632525712 0.2 3.29632525712 # NH3 + 12 0.12 3.02905564168 0.12 2.49451641079 # O + 13 0.12 3.02905564168 0.12 3.02905564168 # OC + +Atoms + + 1 1 11 -0.3 0.0088076654 -0.0395361015 -0.0125765907 # NH3 + 2 1 2 0.33 -0.3781208354 -1.0038773849 -0.01724272 # HC + 3 1 2 0.33 -0.3448285543 0.4901827566 -0.8403800387 # HC + 4 1 2 0.33 -0.3306420078 0.4732826156 0.8294424358 # HC + 5 1 7 0.13 1.526230489 -0.0164860529 -0.0402820599 # CT2 + 6 1 4 0.09 1.8596639218 -0.5263587482 -0.9333647137 # HB + 7 1 4 0.09 1.8904342902 -0.4510777655 0.8809945249 # HB + 8 1 5 0.51 2.0135471936 1.4020233344 -0.1137107587 # C + 9 1 12 -0.51 1.1818164992 2.2781068718 -0.313197467 # O + 10 2 10 -0.47 3.3194424268 1.666672014 0.0713249543 # NH1 + 11 2 1 0.31 4.0409429153 0.9783582555 0.1693897169 # H + 12 2 6 0.07 3.8529467728 3.0161771739 0.0690818527 # CT1 + 13 2 4 0.09 3.5315723829 3.5134434764 -0.8378950509 # HB + 14 2 8 -0.27 3.3981217437 3.8178727883 1.3071600161 # CT3 + 15 2 3 0.09 2.2921877645 3.9163330652 1.3111957959 # HA + 16 2 3 0.09 3.7142964083 3.3104118688 2.2424399743 # HA + 17 2 3 0.09 3.8238672849 4.8446179852 1.2978564897 # HA + 18 2 5 0.51 5.3731873136 2.946062408 0.0042414688 # C + 19 2 12 -0.51 5.9268816775 1.8488403755 0.0011766847 # O + 20 3 10 -0.47 6.0445386951 4.1111769168 -0.0636349698 # NH1 + 21 3 1 0.31 5.5768514013 4.9972489389 -0.0653583036 # H + 22 3 7 -0.02 7.4911642129 4.2772062416 -0.0346088891 # CT2 + 23 3 4 0.09 7.9218377956 3.8786283786 -0.9407244447 # HB + 24 3 4 0.09 7.8877470687 3.8607990619 0.8803359097 # HB + 25 3 5 0.51 7.7351287723 5.7597169026 -0.0088736611 # C + 26 3 12 -0.51 6.7657462979 6.5189223559 0.005150418 # O + 27 4 10 -0.47 9.0014031015 6.1940849758 -0.0045995102 # NH1 + 28 4 1 0.31 9.8272711181 5.6225413025 -0.008448093 # H + 29 4 7 -0.02 9.4067557288 7.5845263022 0.0016383819 # CT2 + 30 4 4 0.09 9.0736276253 8.0578903151 0.9119901724 # HB + 31 4 4 0.09 9.0736165596 8.064051906 -0.9055419386 # HB + 32 4 9 0.34 10.9382207556 7.612479283 0.0008762597 # CC + 33 4 13 -0.67 11.5487033003 6.5062943609 -0.0007524693 # OC + 34 4 13 -0.67 11.5055524841 8.734223435 0.0013684148 # OC + +Bond Coeffs + + 1 250 1.49 # C CT1 + 2 250 1.49 # C CT2 + 3 370 1.345 # C NH1 + 4 620 1.23 # C O + 5 200 1.522 # CC CT2 + 6 525 1.26 # CC OC + 7 222.5 1.538 # CT1 CT3 + 8 330 1.08 # CT1 HB + 9 320 1.43 # CT1 NH1 + 10 330 1.08 # CT2 HB + 11 320 1.43 # CT2 NH1 + 12 200 1.48 # CT2 NH3 + 13 322 1.111 # CT3 HA + 14 440 0.997 # H NH1 + 15 403 1.04 # HC NH3 + +Bonds + + 1 15 2 1 # HC NH3 + 2 15 3 1 # HC NH3 + 3 15 4 1 # HC NH3 + 4 12 1 5 # CT2 NH3 + 5 2 8 5 # C CT2 + 6 3 8 10 # C NH1 + 7 10 5 6 # CT2 HB + 8 10 5 7 # CT2 HB + 9 4 9 8 # C O + 10 7 14 12 # CT1 CT3 + 11 14 10 11 # H NH1 + 12 9 10 12 # CT1 NH1 + 13 1 18 12 # C CT1 + 14 3 18 20 # C NH1 + 15 8 12 13 # CT1 HB + 16 13 14 15 # CT3 HA + 17 13 14 16 # CT3 HA + 18 13 14 17 # CT3 HA + 19 4 19 18 # C O + 20 14 20 21 # H NH1 + 21 11 20 22 # CT2 NH1 + 22 2 25 22 # C CT2 + 23 3 25 27 # C NH1 + 24 10 22 23 # CT2 HB + 25 10 22 24 # CT2 HB + 26 4 26 25 # C O + 27 14 27 28 # H NH1 + 28 11 27 29 # CT2 NH1 + 29 5 32 29 # CC CT2 + 30 10 29 30 # CT2 HB + 31 10 29 31 # CT2 HB + 32 6 32 34 # CC OC + 33 6 32 33 # CC OC + +Angle Coeffs + + 1 52 108 0 0 # C CT1 CT3 + 2 50 109.5 0 0 # C CT1 HB + 3 50 107 0 0 # C CT1 NH1 + 4 50 109.5 0 0 # C CT2 HB + 5 50 107 0 0 # C CT2 NH1 + 6 43.7 110 0 0 # C CT2 NH3 + 7 50 120 0 0 # C NH1 CT1 + 8 50 120 0 0 # C NH1 CT2 + 9 34 123 0 0 # C NH1 H + 10 50 109.5 0 0 # CC CT2 HB + 11 50 107 0 0 # CC CT2 NH1 + 12 80 116.5 0 0 # CT1 C NH1 + 13 80 121 0 0 # CT1 C O + 14 33.43 110.1 22.53 2.179 # CT1 CT3 HA + 15 35 117 0 0 # CT1 NH1 H + 16 80 116.5 0 0 # CT2 C NH1 + 17 80 121 0 0 # CT2 C O + 18 40 118 50 2.388 # CT2 CC OC + 19 35 117 0 0 # CT2 NH1 H + 20 30 109.5 20 2.074 # CT2 NH3 HC + 21 35 111 0 0 # CT3 CT1 HB + 22 70 113.5 0 0 # CT3 CT1 NH1 + 23 35.5 108.4 5.4 1.802 # HA CT3 HA + 24 48 108 0 0 # HB CT1 NH1 + 25 36 115 0 0 # HB CT2 HB + 26 48 108 0 0 # HB CT2 NH1 + 27 51.5 107.5 0 0 # HB CT2 NH3 + 28 44 109.5 0 0 # HC NH3 HC + 29 80 122.5 0 0 # NH1 C O + 30 100 124 70 2.225 # OC CC OC + +Angles + + 1 28 2 1 3 # HC NH3 HC + 2 28 2 1 4 # HC NH3 HC + 3 20 2 1 5 # CT2 NH3 HC + 4 28 3 1 4 # HC NH3 HC + 5 20 3 1 5 # CT2 NH3 HC + 6 20 4 1 5 # CT2 NH3 HC + 7 27 1 5 6 # HB CT2 NH3 + 8 27 1 5 7 # HB CT2 NH3 + 9 6 1 5 8 # C CT2 NH3 + 10 25 6 5 7 # HB CT2 HB + 11 4 6 5 8 # C CT2 HB + 12 4 7 5 8 # C CT2 HB + 13 17 5 8 9 # CT2 C O + 14 16 5 8 10 # CT2 C NH1 + 15 29 9 8 10 # NH1 C O + 16 9 8 10 11 # C NH1 H + 17 7 8 10 12 # C NH1 CT1 + 18 15 11 10 12 # CT1 NH1 H + 19 24 10 12 13 # HB CT1 NH1 + 20 22 10 12 14 # CT3 CT1 NH1 + 21 3 10 12 18 # C CT1 NH1 + 22 21 13 12 14 # CT3 CT1 HB + 23 2 13 12 18 # C CT1 HB + 24 1 14 12 18 # C CT1 CT3 + 25 14 12 14 15 # CT1 CT3 HA + 26 14 12 14 16 # CT1 CT3 HA + 27 14 12 14 17 # CT1 CT3 HA + 28 23 15 14 16 # HA CT3 HA + 29 23 15 14 17 # HA CT3 HA + 30 23 16 14 17 # HA CT3 HA + 31 13 12 18 19 # CT1 C O + 32 12 12 18 20 # CT1 C NH1 + 33 29 19 18 20 # NH1 C O + 34 9 18 20 21 # C NH1 H + 35 8 18 20 22 # C NH1 CT2 + 36 19 21 20 22 # CT2 NH1 H + 37 26 20 22 23 # HB CT2 NH1 + 38 26 20 22 24 # HB CT2 NH1 + 39 5 20 22 25 # C CT2 NH1 + 40 25 23 22 24 # HB CT2 HB + 41 4 23 22 25 # C CT2 HB + 42 4 24 22 25 # C CT2 HB + 43 17 22 25 26 # CT2 C O + 44 16 22 25 27 # CT2 C NH1 + 45 29 26 25 27 # NH1 C O + 46 9 25 27 28 # C NH1 H + 47 8 25 27 29 # C NH1 CT2 + 48 19 28 27 29 # CT2 NH1 H + 49 26 27 29 30 # HB CT2 NH1 + 50 26 27 29 31 # HB CT2 NH1 + 51 11 27 29 32 # CC CT2 NH1 + 52 25 30 29 31 # HB CT2 HB + 53 10 30 29 32 # CC CT2 HB + 54 10 31 29 32 # CC CT2 HB + 55 18 29 32 33 # CT2 CC OC + 56 18 29 32 34 # CT2 CC OC + 57 30 33 32 34 # OC CC OC + +Dihedral Coeffs + + 1 0.2 3 0 1 # C CT1 CT3 HA + 2 0.2 1 180 1 # C CT1 NH1 C + 3 0 1 0 1 # C CT1 NH1 H + 4 0.2 1 180 1 # C CT2 NH1 C + 5 0 1 0 1 # C CT2 NH1 H + 6 0.1 3 0 1 # C CT2 NH3 HC + 7 1.8 1 0 1 # C NH1 CT1 CT3 + 8 0 1 0 1 # C NH1 CT1 HB + 9 0.2 1 180 1 # C NH1 CT2 CC + 10 0 1 0 1 # C NH1 CT2 HB + 11 0 1 0 1 # CC CT2 NH1 H + 12 1.6 1 0 1 # CT1 C NH1 CT2 + 13 2.5 2 180 0 # CT1 C NH1 CT2 + 14 2.5 2 180 1 # CT1 C NH1 H + 15 1.6 1 0 1 # CT1 NH1 C CT2 + 16 2.5 2 180 0 # CT1 NH1 C CT2 + 17 2.5 2 180 1 # CT1 NH1 C O + 18 1.6 1 0 1 # CT2 C NH1 CT2 + 19 2.5 2 180 0 # CT2 C NH1 CT2 + 20 2.5 2 180 1 # CT2 C NH1 H + 21 2.5 2 180 1 # CT2 NH1 C O + 22 0 1 0 1 # CT3 CT1 C NH1 + 23 1.4 1 0 1 # CT3 CT1 C O + 24 0 1 0 1 # CT3 CT1 NH1 H + 25 2.5 2 180 1 # H NH1 C O + 26 0 1 0 1 # H NH1 CT1 HB + 27 0 1 0 1 # H NH1 CT2 HB + 28 0.2 3 0 1 # HA CT3 CT1 HB + 29 0.2 3 0 1 # HA CT3 CT1 NH1 + 30 0 1 0 1 # HB CT1 C NH1 + 31 0 1 0 1 # HB CT1 C O + 32 0 1 0 1 # HB CT2 C NH1 + 33 0 1 0 1 # HB CT2 C O + 34 0.05 6 180 1 # HB CT2 CC OC + 35 0.1 3 0 1 # HB CT2 NH3 HC + 36 0.6 1 0 1 # NH1 C CT1 NH1 + 37 0.6 1 0 1 # NH1 C CT2 NH1 + 38 0.4 1 0 1 # NH1 C CT2 NH3 + 39 0 1 0 1 # NH1 CT1 C O + 40 0 1 0 1 # NH1 CT2 C O + 41 0.05 6 180 1 # NH1 CT2 CC OC + 42 0 1 0 1 # NH3 CT2 C O + +Dihedrals + + 1 42 1 5 8 9 # NH3 CT2 C O + 2 38 1 5 8 10 # NH1 C CT2 NH3 + 3 35 2 1 5 6 # HB CT2 NH3 HC + 4 35 2 1 5 7 # HB CT2 NH3 HC + 5 6 2 1 5 8 # C CT2 NH3 HC + 6 35 3 1 5 6 # HB CT2 NH3 HC + 7 35 3 1 5 7 # HB CT2 NH3 HC + 8 6 3 1 5 8 # C CT2 NH3 HC + 9 35 4 1 5 6 # HB CT2 NH3 HC + 10 35 4 1 5 7 # HB CT2 NH3 HC + 11 6 4 1 5 8 # C CT2 NH3 HC + 12 20 5 8 10 11 # CT2 C NH1 H + 13 15 5 8 10 12 # CT1 NH1 C CT2 + 14 16 5 8 10 12 # CT1 NH1 C CT2 + 15 33 6 5 8 9 # HB CT2 C O + 16 32 6 5 8 10 # HB CT2 C NH1 + 17 33 7 5 8 9 # HB CT2 C O + 18 32 7 5 8 10 # HB CT2 C NH1 + 19 8 8 10 12 13 # C NH1 CT1 HB + 20 7 8 10 12 14 # C NH1 CT1 CT3 + 21 2 8 10 12 18 # C CT1 NH1 C + 22 25 9 8 10 11 # H NH1 C O + 23 17 9 8 10 12 # CT1 NH1 C O + 24 29 10 12 14 15 # HA CT3 CT1 NH1 + 25 29 10 12 14 16 # HA CT3 CT1 NH1 + 26 29 10 12 14 17 # HA CT3 CT1 NH1 + 27 39 10 12 18 19 # NH1 CT1 C O + 28 36 10 12 18 20 # NH1 C CT1 NH1 + 29 26 11 10 12 13 # H NH1 CT1 HB + 30 24 11 10 12 14 # CT3 CT1 NH1 H + 31 3 11 10 12 18 # C CT1 NH1 H + 32 14 12 18 20 21 # CT1 C NH1 H + 33 12 12 18 20 22 # CT1 C NH1 CT2 + 34 13 12 18 20 22 # CT1 C NH1 CT2 + 35 28 13 12 14 15 # HA CT3 CT1 HB + 36 28 13 12 14 16 # HA CT3 CT1 HB + 37 28 13 12 14 17 # HA CT3 CT1 HB + 38 31 13 12 18 19 # HB CT1 C O + 39 30 13 12 18 20 # HB CT1 C NH1 + 40 23 14 12 18 19 # CT3 CT1 C O + 41 22 14 12 18 20 # CT3 CT1 C NH1 + 42 1 15 14 12 18 # C CT1 CT3 HA + 43 1 16 14 12 18 # C CT1 CT3 HA + 44 1 17 14 12 18 # C CT1 CT3 HA + 45 10 18 20 22 23 # C NH1 CT2 HB + 46 10 18 20 22 24 # C NH1 CT2 HB + 47 4 18 20 22 25 # C CT2 NH1 C + 48 25 19 18 20 21 # H NH1 C O + 49 21 19 18 20 22 # CT2 NH1 C O + 50 40 20 22 25 26 # NH1 CT2 C O + 51 37 20 22 25 27 # NH1 C CT2 NH1 + 52 27 21 20 22 23 # H NH1 CT2 HB + 53 27 21 20 22 24 # H NH1 CT2 HB + 54 5 21 20 22 25 # C CT2 NH1 H + 55 20 22 25 27 28 # CT2 C NH1 H + 56 18 22 25 27 29 # CT2 C NH1 CT2 + 57 19 22 25 27 29 # CT2 C NH1 CT2 + 58 33 23 22 25 26 # HB CT2 C O + 59 32 23 22 25 27 # HB CT2 C NH1 + 60 33 24 22 25 26 # HB CT2 C O + 61 32 24 22 25 27 # HB CT2 C NH1 + 62 10 25 27 29 30 # C NH1 CT2 HB + 63 10 25 27 29 31 # C NH1 CT2 HB + 64 9 25 27 29 32 # C NH1 CT2 CC + 65 25 26 25 27 28 # H NH1 C O + 66 21 26 25 27 29 # CT2 NH1 C O + 67 41 27 29 32 33 # NH1 CT2 CC OC + 68 41 27 29 32 34 # NH1 CT2 CC OC + 69 27 28 27 29 30 # H NH1 CT2 HB + 70 27 28 27 29 31 # H NH1 CT2 HB + 71 11 28 27 29 32 # CC CT2 NH1 H + 72 34 30 29 32 33 # HB CT2 CC OC + 73 34 30 29 32 34 # HB CT2 CC OC + 74 34 31 29 32 33 # HB CT2 CC OC + 75 34 31 29 32 34 # HB CT2 CC OC + +Improper Coeffs + + 1 120 0 # C CT1 NH1 O + 2 120 0 # C CT2 NH1 O + 3 96 0 # CC CT2 OC OC + 4 20 0 # H CT1 C NH1 + 5 20 0 # H CT2 C NH1 + +Impropers + + 1 2 8 5 10 9 # C CT2 NH1 O + 2 4 10 8 12 11 # H CT1 C NH1 + 3 1 18 12 20 19 # C CT1 NH1 O + 4 5 20 18 22 21 # H CT2 C NH1 + 5 2 25 22 27 26 # C CT2 NH1 O + 6 5 27 25 29 28 # H CT2 C NH1 + 7 3 32 29 34 33 # CC CT2 OC OC + +CMAP + + 1 1 8 10 12 18 20 + 2 5 18 20 22 25 27 diff --git a/examples/cmap/in.cmap b/examples/cmap/in.cmap new file mode 100644 index 0000000000..d2b2714b82 --- /dev/null +++ b/examples/cmap/in.cmap @@ -0,0 +1,36 @@ +# Created by charmm2lammps v1.8.2.6 beta on Thu Mar 3 20:56:57 EST 2016 + +units real +neigh_modify delay 2 every 1 +#newton off + +boundary p p p + +atom_style full +bond_style harmonic +angle_style charmm +dihedral_style charmm +improper_style harmonic + +pair_style lj/charmm/coul/charmm 8 12 +#pair_style lj/charmmfsw/coul/charmmfsh 8 12 +pair_modify mix arithmetic + +fix cmap all cmap charmm22.cmap +fix_modify cmap energy yes + +read_data gagg.data fix cmap crossterm CMAP + +special_bonds charmm +fix 1 all nve + +#fix 1 all nvt temp 300 300 100.0 +#fix 2 all shake 1e-9 500 0 m 1.0 + +velocity all create 0.0 12345678 dist uniform + +thermo 1000 +thermo_style custom step ecoul evdwl ebond eangle edihed f_cmap eimp +timestep 2.0 + +run 100000 diff --git a/examples/cmap/log.23Sep16.cmap.g++.1 b/examples/cmap/log.23Sep16.cmap.g++.1 new file mode 100644 index 0000000000..fe6caed913 --- /dev/null +++ b/examples/cmap/log.23Sep16.cmap.g++.1 @@ -0,0 +1,200 @@ +LAMMPS (21 Sep 2016) +# Created by charmm2lammps v1.8.2.6 beta on Thu Mar 3 20:56:57 EST 2016 + +units real +neigh_modify delay 2 every 1 +#newton off + +boundary p p p + +atom_style full +bond_style harmonic +angle_style charmm +dihedral_style charmm +improper_style harmonic + +pair_style lj/charmm/coul/charmm 8 12 +#pair_style lj/charmmfsw/coul/charmmfsh 8 12 +pair_modify mix arithmetic + +fix cmap all cmap charmm22.cmap +fix_modify cmap energy yes + +read_data gagg.data fix cmap crossterm CMAP + orthogonal box = (-34.4147 -36.1348 -39.3491) to (45.5853 43.8652 40.6509) + 1 by 1 by 1 MPI processor grid + reading atoms ... + 34 atoms + scanning bonds ... + 4 = max bonds/atom + scanning angles ... + 6 = max angles/atom + scanning dihedrals ... + 12 = max dihedrals/atom + scanning impropers ... + 1 = max impropers/atom + reading bonds ... + 33 bonds + reading angles ... + 57 angles + reading dihedrals ... + 75 dihedrals + reading impropers ... + 7 impropers + 4 = max # of 1-2 neighbors + 7 = max # of 1-3 neighbors + 13 = max # of 1-4 neighbors + 16 = max # of special neighbors + +special_bonds charmm +fix 1 all nve + +#fix 1 all nvt temp 300 300 100.0 +#fix 2 all shake 1e-9 500 0 m 1.0 + +velocity all create 0.0 12345678 dist uniform + +thermo 1000 +thermo_style custom step ecoul evdwl ebond eangle edihed f_cmap eimp +timestep 2.0 + +run 100000 +Neighbor list info ... + 1 neighbor list requests + update every 1 steps, delay 2 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 14 + ghost atom cutoff = 14 + binsize = 7 -> bins = 12 12 12 +Memory usage per processor = 14.6355 Mbytes +Step E_coul E_vdwl E_bond E_angle E_dihed f_cmap E_impro + 0 26.542777 -0.93822087 1.2470497 4.8441789 4.5432816 -1.473352 0.10453023 + 1000 28.673005 -0.47724367 0.80029132 3.151679 4.4684446 -2.3928648 0.18604953 + 2000 27.67955 -1.170342 0.72018905 4.0400131 4.4713764 -2.5490207 0.21834436 + 3000 29.256656 -0.35856055 0.73303546 3.7411606 4.4710568 -2.8939692 0.37728884 + 4000 30.097549 -1.1353905 0.79007053 3.0688444 4.4091469 -2.3383587 0.20743631 + 5000 28.357525 -1.0723742 0.9180297 3.6579424 4.8792664 -2.3185572 0.088366962 + 6000 29.214175 -0.95299225 0.81926009 3.6805429 4.6742897 -2.9343577 0.26697813 + 7000 27.018614 -0.52423475 0.72502764 3.8840137 4.7780956 -2.3916009 0.24952584 + 8000 29.682167 -1.0939711 0.76111486 3.1090116 4.9359719 -2.5662984 0.1411154 + 9000 27.909695 -0.80905986 0.78952533 4.203187 4.1301204 -2.000402 0.088859259 + 10000 27.480298 -0.86273377 1.1293962 4.3857421 4.899282 -3.3895621 0.12126215 + 11000 28.303203 -1.0221152 0.62762348 4.055414 4.5863024 -2.5842816 0.17996907 + 12000 28.311127 -0.94227367 0.91859012 3.6673926 4.7018632 -3.902715 0.30065704 + 13000 30.818607 -1.5220116 0.95710386 3.3364371 4.543427 -3.0423067 0.16712905 + 14000 27.643736 -1.0144117 0.95806952 4.1046912 4.800236 -4.0534389 0.29293405 + 15000 27.660491 -1.0390086 0.78061056 4.1139174 4.7197202 -2.3670379 0.22126985 + 16000 27.845157 -0.63654502 0.78007478 3.9365994 4.949418 -3.1470214 0.22335355 + 17000 28.44772 -1.0255112 0.70402007 4.0573343 4.2887527 -2.2099596 0.048050962 + 18000 27.128323 -0.96218536 1.1327159 4.3222585 4.326607 -2.2881766 0.13491257 + 19000 27.337633 -0.78999574 0.80152298 4.2239689 4.7073478 -2.2924164 0.12710292 + 20000 27.780537 -0.46458072 0.79707671 3.7232618 4.943417 -2.5290628 0.26191223 + 21000 26.435484 -0.7803224 1.0753217 4.4196051 5.9945933 -2.3340925 0.16448475 + 22000 28.619429 -1.1623149 0.9401731 3.8508844 5.1636737 -2.5551846 0.25318434 + 23000 28.399338 -0.79700962 0.85575503 4.488526 4.5975422 -2.5663549 0.13601693 + 24000 29.645532 -1.158744 0.83180313 3.8193399 4.60319 -2.6847864 0.24260466 + 25000 28.695339 -1.4802204 0.76583757 3.6786272 4.8959496 -2.3627896 0.080867326 + 26000 28.149711 -1.029689 0.79383806 3.7885067 4.3345813 -2.1041553 0.14598209 + 27000 29.580373 -1.0525813 1.0262723 3.7767318 4.6119758 -2.2802386 0.088556038 + 28000 28.44308 -0.93411225 0.8794395 3.948079 4.780246 -2.1814583 0.14340149 + 29000 29.335621 -1.6087988 0.71803091 3.7819186 4.6688385 -2.4282242 0.16061111 + 30000 28.706138 -1.3938241 0.67713818 4.031275 4.4756505 -2.1807056 0.11461356 + 31000 27.451944 0.010297225 0.65064883 3.6402029 4.3607811 -2.5511516 0.12637237 + 32000 27.070878 -1.103158 1.1932199 5.1329709 4.5201653 -2.2224479 0.11215427 + 33000 29.889976 -1.6228316 0.69407996 3.5361991 4.3502767 -1.9847454 0.09089949 + 34000 28.223151 -0.927208 1.043253 3.4650939 5.1028142 -2.8127219 0.10648823 + 35000 27.985986 -0.48153861 0.63878449 3.3724641 4.9551679 -2.6565919 0.12123115 + 36000 28.580688 -1.4500694 1.055762 4.0490064 4.423782 -2.3103578 0.072747638 + 37000 29.192947 -0.49678176 1.1146731 2.9233947 4.5738603 -2.4376144 0.22874047 + 38000 26.954594 -0.53812359 0.79230685 4.3356989 5.0284656 -2.3791255 0.0486081 + 39000 27.567555 -0.57870028 0.73614374 4.191991 4.9209556 -2.6122044 0.08635571 + 40000 28.494172 -0.79057135 0.79072816 4.1893209 4.4826919 -2.4179635 0.14612898 + 41000 28.44904 -1.1002948 0.93405654 4.3586358 4.4338415 -2.2950944 0.15705834 + 42000 28.95725 -1.0297067 1.1632348 4.274711 4.9979487 -2.7611464 0.15944725 + 43000 28.640394 -0.70938882 0.68100893 3.1844315 5.1817454 -2.2837487 0.14189233 + 44000 27.997558 -1.0115099 0.59125208 4.0883422 4.6033907 -2.2775964 0.094273258 + 45000 27.67163 -0.67992747 1.1225552 3.9020703 4.8171056 -2.1952679 0.041418433 + 46000 28.822607 -0.6687953 0.74160059 3.3193715 4.5546965 -2.3024572 0.047569065 + 47000 29.20147 -1.4456785 0.79223353 3.8288813 4.5811826 -2.5154936 0.061230141 + 48000 27.843026 -1.0222301 0.87322137 4.3432743 4.4266307 -2.1414153 0.06802794 + 49000 28.199573 -1.1887794 1.2781088 4.0779644 4.5881353 -2.319775 0.094803547 + 50000 28.759212 -1.354416 0.68534569 3.8394841 4.2308134 -2.1281844 0.1395951 + 51000 27.876455 -1.5705462 0.76557156 4.5335223 4.523708 -2.203702 0.14679803 + 52000 27.930587 -1.2277489 0.96071516 3.960953 5.1152188 -2.4101451 0.060949521 + 53000 27.031236 -1.4746477 1.2341141 5.0540975 4.3656865 -2.1288513 0.092725656 + 54000 28.809394 -1.1162427 0.94350207 3.4013958 4.4755547 -2.3342811 0.18855912 + 55000 28.948415 -1.1560418 0.6260139 3.5386373 4.5244978 -2.340212 0.17474657 + 56000 28.048368 -0.95784532 0.76432571 4.1404665 4.4570033 -2.0899628 0.045693628 + 57000 28.707642 -1.366574 0.9907873 3.729903 4.3131997 -2.2777698 0.065420213 + 58000 26.361663 -1.0424403 1.0452563 5.0977108 4.7035231 -2.3101244 0.13671642 + 59000 29.218218 -1.2210564 0.62435875 3.4236327 4.5481681 -2.1575943 0.037984042 + 60000 27.655546 -1.1053224 0.86323501 3.7641375 4.8946898 -2.2422249 0.077725979 + 61000 27.252108 -1.3744824 1.1150806 5.0444848 4.4878135 -2.2743829 0.058331257 + 62000 27.163469 -1.1715781 0.72099321 4.5295501 4.9509918 -2.2993961 0.050401105 + 63000 29.581575 -1.2238537 0.86303245 3.1194038 5.2218965 -2.5002427 0.055032632 + 64000 27.897822 -1.1011516 0.74540883 4.2869228 4.3394269 -2.2552393 0.1403321 + 65000 27.083245 -1.0633392 0.92771724 5.0805224 4.2747962 -2.2388039 0.064196692 + 66000 29.072723 -1.5514209 0.89798805 4.2600224 4.4261812 -2.3524752 0.15067414 + 67000 27.308181 -0.72224802 0.97109517 4.5074578 4.4559352 -2.1381121 0.089297603 + 68000 27.505686 -0.43855431 0.80785812 4.1917251 5.0157721 -2.3382145 0.11105164 + 69000 29.041681 -0.64735378 0.89874684 3.3891579 4.3753361 -2.2320941 0.14716747 + 70000 29.735756 -1.7061457 0.9206878 3.5767878 4.3851664 -2.2516304 0.097196062 + 71000 28.224352 -0.92217702 0.86093586 3.9507157 4.5596589 -2.2173397 0.089116669 + 72000 29.282336 -1.056142 0.65185725 3.8735742 4.4839333 -2.4314756 0.071909704 + 73000 26.257283 -0.64273826 0.98300685 5.063943 5.045958 -2.5544375 0.2180275 + 74000 28.825119 -0.97736616 0.87201848 3.55875 4.3653309 -2.2303567 0.098963875 + 75000 29.239507 -0.96508809 0.74517323 3.4306236 4.7651921 -2.6077732 0.17883654 + 76000 27.349841 -0.50990238 1.1183613 4.4252451 4.4097775 -2.4125794 0.18483606 + 77000 28.130197 -1.4081219 0.94921357 4.2572132 4.5162849 -2.4013797 0.073744606 + 78000 28.235774 -0.9214321 0.6324981 3.8697686 4.8092154 -2.2272847 0.092108346 + 79000 26.732846 -0.55949486 1.0989617 5.0088609 4.4930687 -2.277945 0.03855146 + 80000 28.529208 -0.94244671 0.79407482 3.961106 4.3930011 -2.3127726 0.091124948 + 81000 29.603852 -1.6116062 1.060847 3.7824932 4.151001 -1.9139868 0.19875986 + 82000 28.232876 -1.1833011 1.0182713 3.4195758 5.1394333 -2.4632697 0.28501012 + 83000 29.565482 -1.3479552 0.99056973 3.7851802 4.4781011 -2.7872481 0.2031991 + 84000 28.780274 -1.3073882 1.0512637 4.004638 4.502282 -2.3789146 0.015656202 + 85000 27.262312 -1.1305346 1.203524 4.7938623 4.1747105 -2.0952844 0.054240361 + 86000 28.157348 -1.0662817 0.81163796 3.9912709 4.8320213 -2.255237 0.14698333 + 87000 28.445543 -1.3365026 0.78156195 4.4767689 4.4457575 -2.5008786 0.13879386 + 88000 27.656717 -1.1490599 0.87974869 4.4629952 4.7023033 -2.3258145 0.081904139 + 89000 28.838821 -1.020709 0.85587929 3.7110705 4.4938307 -2.4914483 0.11447952 + 90000 27.356497 -0.59107077 0.81879666 4.5209332 4.4703836 -2.3806717 0.071307775 + 91000 27.780445 -0.80564513 0.94752313 3.8468943 4.2924253 -2.1011134 0.1118672 + 92000 28.555276 -1.3514732 0.80826674 3.9590742 4.5775954 -2.4891232 0.054254978 + 93000 28.747267 -1.2133243 0.75507246 4.1319789 4.9048611 -2.4913887 0.13045693 + 94000 27.479343 -0.69973695 0.99696121 3.5966229 4.549025 -2.4155312 0.41745762 + 95000 27.726945 -1.1905026 1.1120842 4.7433275 4.5386861 -2.7947142 0.33671682 + 96000 28.021114 -1.0341645 0.6663033 4.2397505 4.6203984 -1.9904034 0.10972565 + 97000 28.382022 -1.3916008 1.180588 4.0729621 4.6741792 -2.554927 0.13462346 + 98000 27.895969 -0.7496449 1.3072185 4.2611888 4.3726077 -2.1320701 0.15376665 + 99000 28.517889 -1.2183957 1.279778 3.957647 4.2638434 -2.2888407 0.042705003 + 100000 28.109211 -1.2538948 0.83671785 4.3734766 4.544545 -2.3076497 0.042189096 +Loop time of 2.96683 on 1 procs for 100000 steps with 34 atoms + +Performance: 5824.390 ns/day, 0.004 hours/ns, 33705.963 timesteps/s +100.0% CPU use with 1 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0.98759 | 0.98759 | 0.98759 | 0.0 | 33.29 +Bond | 1.6463 | 1.6463 | 1.6463 | 0.0 | 55.49 +Neigh | 0.007688 | 0.007688 | 0.007688 | 0.0 | 0.26 +Comm | 0.012214 | 0.012214 | 0.012214 | 0.0 | 0.41 +Output | 0.0010295 | 0.0010295 | 0.0010295 | 0.0 | 0.03 +Modify | 0.25684 | 0.25684 | 0.25684 | 0.0 | 8.66 +Other | | 0.05519 | | | 1.86 + +Nlocal: 34 ave 34 max 34 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 0 ave 0 max 0 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 395 ave 395 max 395 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 395 +Ave neighs/atom = 11.6176 +Ave special neighs/atom = 9.52941 +Neighbor list builds = 237 +Dangerous builds = 0 +Total wall time: 0:00:02 diff --git a/examples/cmap/log.23Sep16.cmap.g++.4 b/examples/cmap/log.23Sep16.cmap.g++.4 new file mode 100644 index 0000000000..2d3efa66c6 --- /dev/null +++ b/examples/cmap/log.23Sep16.cmap.g++.4 @@ -0,0 +1,200 @@ +LAMMPS (21 Sep 2016) +# Created by charmm2lammps v1.8.2.6 beta on Thu Mar 3 20:56:57 EST 2016 + +units real +neigh_modify delay 2 every 1 +#newton off + +boundary p p p + +atom_style full +bond_style harmonic +angle_style charmm +dihedral_style charmm +improper_style harmonic + +pair_style lj/charmm/coul/charmm 8 12 +#pair_style lj/charmmfsw/coul/charmmfsh 8 12 +pair_modify mix arithmetic + +fix cmap all cmap charmm22.cmap +fix_modify cmap energy yes + +read_data gagg.data fix cmap crossterm CMAP + orthogonal box = (-34.4147 -36.1348 -39.3491) to (45.5853 43.8652 40.6509) + 1 by 2 by 2 MPI processor grid + reading atoms ... + 34 atoms + scanning bonds ... + 4 = max bonds/atom + scanning angles ... + 6 = max angles/atom + scanning dihedrals ... + 12 = max dihedrals/atom + scanning impropers ... + 1 = max impropers/atom + reading bonds ... + 33 bonds + reading angles ... + 57 angles + reading dihedrals ... + 75 dihedrals + reading impropers ... + 7 impropers + 4 = max # of 1-2 neighbors + 7 = max # of 1-3 neighbors + 13 = max # of 1-4 neighbors + 16 = max # of special neighbors + +special_bonds charmm +fix 1 all nve + +#fix 1 all nvt temp 300 300 100.0 +#fix 2 all shake 1e-9 500 0 m 1.0 + +velocity all create 0.0 12345678 dist uniform + +thermo 1000 +thermo_style custom step ecoul evdwl ebond eangle edihed f_cmap eimp +timestep 2.0 + +run 100000 +Neighbor list info ... + 1 neighbor list requests + update every 1 steps, delay 2 steps, check yes + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 14 + ghost atom cutoff = 14 + binsize = 7 -> bins = 12 12 12 +Memory usage per processor = 15.9307 Mbytes +Step E_coul E_vdwl E_bond E_angle E_dihed f_cmap E_impro + 0 26.542777 -0.93822087 1.2470497 4.8441789 4.5432816 -1.473352 0.10453023 + 1000 28.673005 -0.47724367 0.80029132 3.151679 4.4684446 -2.3928648 0.18604953 + 2000 27.67955 -1.170342 0.72018905 4.0400131 4.4713764 -2.5490207 0.21834436 + 3000 29.256656 -0.35856055 0.73303546 3.7411606 4.4710568 -2.8939692 0.37728884 + 4000 30.097549 -1.1353905 0.79007053 3.0688444 4.4091469 -2.3383587 0.20743631 + 5000 28.357525 -1.0723742 0.9180297 3.6579424 4.8792663 -2.3185572 0.088366962 + 6000 29.214175 -0.95299239 0.81926011 3.6805428 4.6742897 -2.9343578 0.26697816 + 7000 27.018614 -0.52423469 0.72502751 3.8840141 4.7780958 -2.3916014 0.24952572 + 8000 29.682494 -1.0940368 0.76113051 3.1089345 4.9357863 -2.5662256 0.14112613 + 9000 27.853918 -0.7913741 0.79503268 4.2177256 4.146792 -2.00475 0.090585666 + 10000 27.13754 -0.80551128 1.1325023 4.4718283 5.2460631 -3.4947725 0.11893125 + 11000 28.277434 -1.4897448 0.90075953 4.1895717 4.3594269 -1.9553119 0.090222212 + 12000 28.630973 -1.222206 0.67796385 3.3905661 4.9691334 -2.9052721 0.13897658 + 13000 28.593007 -0.95684026 0.75585196 3.7242568 4.7417932 -2.3893117 0.2074121 + 14000 26.147115 -0.6026921 0.93591488 5.1292829 4.9821952 -2.2571835 0.11872421 + 15000 26.29432 -0.82424162 1.048979 4.5569495 5.1189308 -2.9750422 0.16195676 + 16000 29.189992 -0.80998247 0.74093508 3.8299275 4.4536688 -2.5497538 0.19155639 + 17000 25.878012 -0.3519646 1.0988924 4.7359591 5.3923098 -2.7211029 0.13405223 + 18000 27.726135 -0.28229987 0.63072344 4.1777888 4.7237271 -2.2177157 0.15939372 + 19000 27.153504 -0.66477422 0.77910129 4.2036117 5.113851 -2.3494315 0.094793307 + 20000 28.044833 -1.2835827 0.88745367 3.9955526 4.5077788 -3.0116467 0.17197859 + 21000 27.205696 -0.74090037 1.0023251 4.3421733 4.912671 -2.3473271 0.26089356 + 22000 27.385785 -0.93740972 0.84554838 4.562743 4.883866 -2.2110955 0.11573301 + 23000 27.05534 -0.95605442 0.96719024 3.9277618 5.0359014 -2.6135949 0.21368061 + 24000 28.273378 -0.97543103 0.8983443 4.2067985 4.4782971 -2.4230505 0.30311692 + 25000 27.477789 -0.20383849 0.8380706 3.8037992 4.8312504 -2.5831791 0.093843746 + 26000 30.344199 -1.9773473 0.92882437 3.7821405 4.5176677 -2.3020968 0.2194307 + 27000 27.32767 -0.9803839 0.92988865 3.7611603 5.0328211 -2.4647656 0.18213622 + 28000 27.34208 -1.037938 0.74488346 4.1727342 4.7056812 -2.2718346 0.17741362 + 29000 27.682777 -0.51006495 0.57074224 4.7332237 4.7080462 -2.0491512 0.2130517 + 30000 24.925731 0.13670248 0.84976065 4.4143762 6.0677158 -3.5479173 0.28059419 + 31000 28.623419 -0.90725708 1.0710501 3.6930688 4.6639301 -2.2225373 0.20988139 + 32000 27.732286 -1.1948367 0.89230134 4.4398373 4.8923907 -3.5849327 0.49167488 + 33000 28.800772 -1.5319589 0.93455495 4.1634728 4.6107706 -2.3503486 0.22636535 + 34000 27.374398 -1.0957453 0.89450276 3.9829508 4.991786 -2.3548834 0.15869465 + 35000 28.38753 -0.89261166 0.90000776 3.536864 4.4293294 -2.4218118 0.10640557 + 36000 27.713974 0.088038031 0.85190574 3.8969601 4.6256355 -2.7935475 0.34671662 + 37000 29.13007 -1.378597 0.74412556 3.131538 4.6458653 -2.9373734 0.38035616 + 38000 28.556573 -1.4055344 1.139984 4.0035753 4.2938358 -2.489329 0.25338326 + 39000 26.447036 -1.1829705 0.87032438 5.0804461 4.5772023 -2.7346466 0.32165802 + 40000 27.991454 -0.64295679 0.61020872 4.165871 4.4623087 -2.2244194 0.13826991 + 41000 29.483296 -1.2400745 0.66926627 3.3473666 4.5766617 -2.3051145 0.12171554 + 42000 26.948627 -1.2162288 1.1440628 4.3993073 5.1176533 -2.4734485 0.15497709 + 43000 28.04459 -0.26543193 0.83647367 3.5160747 4.6964397 -2.2805068 0.12618821 + 44000 28.213608 -1.216128 0.9132792 4.0206483 4.9483599 -2.3387049 0.10132022 + 45000 28.283506 -1.0390766 0.86113772 4.504509 4.7209088 -2.3043085 0.14588362 + 46000 27.433853 -0.57912107 0.78448334 4.5998579 5.1181394 -2.6165094 0.18722528 + 47000 27.552939 -1.1128925 0.80087638 4.3448001 4.8062869 -2.4296883 0.2702479 + 48000 28.874034 -1.3242519 0.71770727 3.5648565 4.4671824 -2.2608958 0.16115978 + 49000 29.216186 -1.2210307 0.76937497 3.9260628 4.7550577 -2.7316081 0.085505664 + 50000 28.065856 -1.1545547 0.86953819 4.4137666 4.732157 -2.4450867 0.23320539 + 51000 26.308975 -0.99728352 0.90408444 4.2400186 5.6340425 -2.2090554 0.079882158 + 52000 28.517571 -1.5027398 0.83520278 3.8176552 4.3001251 -2.0731682 0.1665375 + 53000 28.77579 -1.3564268 0.97253881 3.6866407 4.8532347 -2.5330776 0.17668411 + 54000 29.135315 -1.0994106 0.67605671 3.6819254 4.3134408 -1.9796929 0.076951331 + 55000 26.168938 -0.76247492 0.88784685 4.6533473 6.0484793 -2.1334561 0.036876985 + 56000 27.471775 -0.68648837 1.0576168 4.0354311 4.4767052 -2.2368959 0.24950568 + 57000 29.787083 -1.4914384 1.0702944 3.5388133 4.5173097 -2.6694464 0.27937092 + 58000 28.705448 -1.3016617 0.63337853 3.9552713 4.4119825 -1.8774657 0.17540021 + 59000 29.130155 -0.91647363 0.84384883 3.1076903 4.5346348 -2.3457338 0.16674486 + 60000 26.874199 -0.81598034 1.3432151 5.1322624 4.9545484 -2.9566615 0.25950486 + 61000 27.401306 -0.82895856 1.1636949 4.020154 4.5745928 -2.601466 0.18061051 + 62000 28.930313 -1.5231967 0.85173243 4.3517328 4.4878662 -2.5859205 0.1755493 + 63000 26.56874 0.026147233 0.60836216 4.4231618 4.4390677 -2.1721849 0.08594237 + 64000 26.729023 -0.76953985 0.76734633 4.5104288 5.0886456 -2.2118551 0.11339216 + 65000 28.900471 -1.3901477 0.86194657 4.2774976 4.498325 -2.3672362 0.20668335 + 66000 26.884253 -0.21198879 0.98509625 4.0843117 4.4344172 -2.3289416 0.23631017 + 67000 27.210888 -0.84075559 1.0396559 4.7253607 4.4314589 -2.2985702 0.19326507 + 68000 28.042102 -1.1898715 1.053534 3.8748712 4.4358449 -2.3998723 0.2431659 + 69000 28.939141 -1.6968936 0.98155912 4.0460838 5.0075204 -2.5547087 0.28645131 + 70000 27.15577 -0.85202797 1.1469079 4.7645212 4.6133209 -2.3410451 0.086576572 + 71000 25.507417 -0.27780727 0.95157881 4.8759406 4.853401 -2.9598705 0.41011008 + 72000 29.804703 -1.4847015 0.96345767 3.6797304 4.3678377 -2.4594626 0.14480206 + 73000 28.602798 -1.4906143 0.72497266 4.2442974 4.5360598 -2.3621638 0.14385651 + 74000 28.4928 -0.91319873 1.0377472 3.8033127 4.3991601 -2.4051911 0.095567428 + 75000 26.38168 -0.70733237 1.1557817 5.697939 4.5935618 -2.4285007 0.058980519 + 76000 27.16626 -0.83631031 0.84844246 4.7460887 4.5801472 -2.1260014 0.12845946 + 77000 29.040661 -1.3089499 0.80285084 4.664804 4.5215895 -2.6861939 0.13215598 + 78000 27.477871 -1.0600977 0.88595045 4.6264017 5.4095605 -2.474411 0.10987174 + 79000 26.151797 -0.55779685 0.91382436 4.99964 4.9184022 -2.2547241 0.22854038 + 80000 28.14523 -0.54460026 0.8982411 3.5374555 4.3785673 -2.3196807 0.088567964 + 81000 29.029941 -1.6467789 0.79042284 3.7269899 4.7407998 -2.3795824 0.1408727 + 82000 27.920287 -0.72798032 1.0076975 3.4449461 4.5621371 -2.8239074 0.25103454 + 83000 29.131054 -1.114367 0.76887285 3.459639 4.5163922 -2.607825 0.19991648 + 84000 28.249768 -0.69944068 1.0510846 4.0436296 4.6430538 -2.4213355 0.077299966 + 85000 28.06888 -0.62132922 0.91829312 4.1294147 4.3099557 -2.354063 0.15866186 + 86000 28.664264 -1.1022906 0.87831695 4.5773522 4.6045802 -2.9206875 0.33950063 + 87000 27.960967 -1.2852756 0.77694253 3.9011301 4.9114139 -3.2374868 0.3068138 + 88000 27.190678 -1.2803268 1.1545301 4.5769709 5.2404761 -2.3825838 0.10356039 + 89000 26.792931 -0.44516641 1.0236244 4.2007253 4.7098685 -2.3608551 0.034447062 + 90000 27.173991 -0.87185611 1.065719 4.1953618 4.6856408 -2.6539232 0.16957757 + 91000 28.626528 -1.239257 0.89524651 4.7048012 4.6344201 -2.7367901 0.43534143 + 92000 27.661812 -1.109044 0.92817391 5.0294489 4.3890711 -2.4108669 0.12570139 + 93000 28.156793 -1.0820907 0.92812693 4.938385 4.4901426 -2.4023366 0.30135781 + 94000 28.842149 -1.3524969 1.1451109 4.3125908 4.6959035 -2.6747199 0.2254607 + 95000 27.862247 -1.2119045 1.0218976 4.2614082 4.4931316 -2.6902934 0.16345201 + 96000 27.084973 -0.93738328 1.3984324 4.5647189 4.4232205 -2.2834097 0.11217888 + 97000 27.587078 -0.89397255 0.78218462 3.8944421 4.3981479 -2.4205318 0.16570942 + 98000 27.981746 -1.2380545 0.84847869 4.311441 4.7340377 -2.4270441 0.023565612 + 99000 27.476625 -0.8569146 0.82550381 4.1656963 4.4064921 -2.4169708 0.160814 + 100000 26.121325 -0.63610855 1.0803389 4.9257118 4.7073263 -2.4010334 0.066303044 +Loop time of 2.60487 on 4 procs for 100000 steps with 34 atoms + +Performance: 6633.735 ns/day, 0.004 hours/ns, 38389.667 timesteps/s +99.1% CPU use with 4 MPI tasks x no OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 0.066848 | 0.26055 | 0.6843 | 48.6 | 10.00 +Bond | 0.067332 | 0.45486 | 0.93545 | 55.1 | 17.46 +Neigh | 0.0078266 | 0.007863 | 0.0078835 | 0.0 | 0.30 +Comm | 0.41829 | 1.3207 | 1.8951 | 50.8 | 50.70 +Output | 0.0033038 | 0.0036355 | 0.0040481 | 0.4 | 0.14 +Modify | 0.040861 | 0.15162 | 0.27091 | 27.3 | 5.82 +Other | | 0.4057 | | | 15.57 + +Nlocal: 8.5 ave 15 max 2 min +Histogram: 1 1 0 0 0 0 0 0 1 1 +Nghost: 25.5 ave 32 max 19 min +Histogram: 1 1 0 0 0 0 0 0 1 1 +Neighs: 98.75 ave 257 max 18 min +Histogram: 1 1 1 0 0 0 0 0 0 1 + +Total # of neighbors = 395 +Ave neighs/atom = 11.6176 +Ave special neighs/atom = 9.52941 +Neighbor list builds = 294 +Dangerous builds = 0 +Total wall time: 0:00:02 diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp new file mode 100644 index 0000000000..39081168d8 --- /dev/null +++ b/src/MOLECULE/fix_cmap.cpp @@ -0,0 +1,1432 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Implementation of the CHARMM CMAP; adds an extra energy term for the + peptide backbone dihedrals. The tools/ch2lmp.pl conversion script, which + generates an extra section in the LAMMPS data file, is needed in order to + generate the info used by this fix style. + + Contributing authors: + Xiaohu Hu, CMB/ORNL (hux2@ornl.gov) + David Hyde-Volpe, Tigran Abramyan, and Robert A. Latour (Clemson University) + Chris Lorenz (Kings College-London) + + References: + - MacKerell et al., J. Am. Chem. Soc. 126(2004):698-699. + - MacKerell et al., J. Comput. Chem. 25(2004):1400-1415. + -------------------------------------------------------------------------*/ + +#include "mpi.h" +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "stdio.h" +#include "fix_cmap.h" +#include "atom.h" +#include "atom_vec.h" +#include "update.h" +#include "respa.h" +#include "modify.h" +#include "domain.h" +#include "force.h" +#include "group.h" +#include "comm.h" +#include "math_const.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace MathConst; + +#define MAXLINE 256 +#define LISTDELTA 10000 +#define LB_FACTOR 1.5 + +#define CMAPMAX 6 // max # of CMAP terms stored by one atom +#define CMAPDIM 24 // grid map dimension is 24 x 24 +#define CMAPXMIN -360.0 +#define CMAPXMIN2 -180.0 +#define CMAPDX 15.0 // 360/CMAPDIM + +/* ---------------------------------------------------------------------- */ + +FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +{ + if (narg != 4) error->all(FLERR,"Illegal fix cmap command"); + + restart_global = 1; + restart_peratom = 1; + peatom_flag = 1; + virial_flag = 1; + peratom_freq = 1; + scalar_flag = 1; + global_freq = 1; + extscalar = 1; + extvector = 1; + wd_header = 1; + wd_section = 1; + + MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); + + // allocate memory for CMAP data + + memory->create(g_axis,CMAPDIM,"cmap:g_axis"); + memory->create(cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:grid"); + memory->create(d1cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d1grid"); + memory->create(d2cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d2grid"); + memory->create(d12cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d12grid"); + + // read, initialize, broadcast cmapgrid + + if (me == 0) read_grid_map(arg[3]); + MPI_Bcast(&cmapgrid[0][0][0],6*CMAPDIM*CMAPDIM,MPI_DOUBLE,0,world); + + // perform initial allocation of atom-based arrays + // register with Atom class + + num_crossterm = NULL; + crossterm_type = NULL; + crossterm_atom1 = NULL; + crossterm_atom2 = NULL; + crossterm_atom3 = NULL; + crossterm_atom4 = NULL; + crossterm_atom5 = NULL; + + nmax_previous = 0; + grow_arrays(atom->nmax); + atom->add_callback(0); + atom->add_callback(1); + + // local list of crossterms + + ncmap = 0; + maxcrossterm = 0; + crosstermlist = NULL; +} + +/* --------------------------------------------------------------------- */ + +FixCMAP::~FixCMAP() +{ + // unregister callbacks to this fix from Atom class + + atom->delete_callback(id,0); + atom->delete_callback(id,1); + + memory->destroy(g_axis); + memory->destroy(cmapgrid); + memory->destroy(d1cmapgrid); + memory->destroy(d2cmapgrid); + memory->destroy(d12cmapgrid); + + memory->destroy(crosstermlist); + + memory->destroy(num_crossterm); + memory->destroy(crossterm_type); + memory->destroy(crossterm_atom1); + memory->destroy(crossterm_atom2); + memory->destroy(crossterm_atom3); + memory->destroy(crossterm_atom4); + memory->destroy(crossterm_atom5); +} + +/* ---------------------------------------------------------------------- */ + +int FixCMAP::setmask() +{ + int mask = 0; + mask |= PRE_NEIGHBOR; + mask |= PRE_REVERSE; + mask |= POST_FORCE; + mask |= THERMO_ENERGY; + mask |= POST_FORCE_RESPA; + mask |= MIN_POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixCMAP::init() +{ + int i; + double angle; + + i = 0; + angle = -180.0; + while (angle < 180.0) { + g_axis[i] = angle; + angle += CMAPDX; + i++; + } + + // pre-compute the derivatives of the maps + + for (i = 0; i < 6; i++) + set_map_derivatives(cmapgrid[i],d1cmapgrid[i],d2cmapgrid[i],d12cmapgrid[i]); + + // define newton_bond here in case restart file was read (not data file) + + newton_bond = force->newton_bond; +} + +/* --------------------------------------------------------------------- */ + +void FixCMAP::setup(int vflag) +{ + pre_neighbor(); + + if (strstr(update->integrate_style,"verlet")) + post_force(vflag); + else { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(vflag,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + } +} + +/* --------------------------------------------------------------------- */ + +void FixCMAP::setup_pre_neighbor() +{ + pre_neighbor(); +} + +/* --------------------------------------------------------------------- */ + +void FixCMAP::min_setup(int vflag) +{ + pre_neighbor(); + post_force(vflag); +} + +/* ---------------------------------------------------------------------- + store local neighbor list as if newton_bond = OFF, even if actually ON +------------------------------------------------------------------------- */ + +void FixCMAP::pre_neighbor() +{ + int i,m,itype,atom1,atom2,atom3,atom4,atom5; + + // guesstimate initial length of local crossterm list + // if ncmap was not set (due to read_restart, no read_data), + // then list will grow by LISTDELTA chunks + + if (maxcrossterm == 0) { + if (nprocs == 1) maxcrossterm = ncmap; + else maxcrossterm = static_cast (LB_FACTOR*ncmap/nprocs); + memory->create(crosstermlist,maxcrossterm,6,"cmap:crosstermlist"); + } + + int nlocal = atom->nlocal; + + ncrosstermlist = 0; + + for (i = 0; i < nlocal; i++) { + for (m = 0; m < num_crossterm[i]; m++) { + atom1 = atom->map(crossterm_atom1[i][m]); + atom2 = atom->map(crossterm_atom2[i][m]); + atom3 = atom->map(crossterm_atom3[i][m]); + atom4 = atom->map(crossterm_atom4[i][m]); + atom5 = atom->map(crossterm_atom5[i][m]); + + if (atom1 == -1 || atom2 == -1 || atom3 == -1 || + atom4 == -1 || atom5 == -1) { + char str[128]; + sprintf(str,"CMAP atoms " + TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " + TAGINT_FORMAT " " TAGINT_FORMAT + " missing on proc %d at step " BIGINT_FORMAT, + crossterm_atom1[i][m],crossterm_atom2[i][m], + crossterm_atom3[i][m],crossterm_atom4[i][m], + crossterm_atom5[i][m],me,update->ntimestep); + error->one(FLERR,str); + } + atom1 = domain->closest_image(i,atom1); + atom2 = domain->closest_image(i,atom2); + atom3 = domain->closest_image(i,atom3); + atom4 = domain->closest_image(i,atom4); + atom5 = domain->closest_image(i,atom5); + + if (i <= atom1 && i <= atom2 && i <= atom3 && + i <= atom4 && i <= atom5) { + if (ncrosstermlist == maxcrossterm) { + maxcrossterm += LISTDELTA; + memory->grow(crosstermlist,maxcrossterm,6,"cmap:crosstermlist"); + } + crosstermlist[ncrosstermlist][0] = atom1; + crosstermlist[ncrosstermlist][1] = atom2; + crosstermlist[ncrosstermlist][2] = atom3; + crosstermlist[ncrosstermlist][3] = atom4; + crosstermlist[ncrosstermlist][4] = atom5; + crosstermlist[ncrosstermlist][5] = crossterm_type[i][m]; + ncrosstermlist++; + } + } + } +} + +/* ---------------------------------------------------------------------- + store eflag, so can use it in post_force to tally per-atom energies +------------------------------------------------------------------------- */ + +void FixCMAP::pre_reverse(int eflag, int vflag) +{ + eflag_caller = eflag; +} + +/* ---------------------------------------------------------------------- + compute CMAP terms as if newton_bond = OFF, even if actually ON +------------------------------------------------------------------------- */ + +void FixCMAP::post_force(int vflag) +{ + int n,i1,i2,i3,i4,i5,type,nlist; + int li1, li2, mli1,mli2,mli11,mli21,t1,li3,li4,mli3,mli4,mli31,mli41; + int list[5]; + // vectors needed to calculate the cross-term dihedral angles + double vb21x,vb21y,vb21z,vb32x,vb32y,vb32z,vb34x,vb34y,vb34z; + double vb23x,vb23y,vb23z; + double vb43x,vb43y,vb43z,vb45x,vb45y,vb45z,a1x,a1y,a1z,b1x,b1y,b1z; + double a2x,a2y,a2z,b2x,b2y,b2z,r32,a1sq,b1sq,a2sq,b2sq,dpr21r32,dpr34r32; + double dpr32r43,dpr45r43,r43,vb12x,vb12y,vb12z,vb54x,vb54y,vb54z; + // cross-term dihedral angles + double phi,psi,phi1,psi1; + double f1[3],f2[3],f3[3],f4[3],f5[3],vcmap[6]; + double gs[4],d1gs[4],d2gs[4],d12gs[4]; + double engfraction; + // vectors needed for the gradient/force calculation + double dphidr1x,dphidr1y,dphidr1z,dphidr2x,dphidr2y,dphidr2z; + double dphidr3x,dphidr3y,dphidr3z,dphidr4x,dphidr4y,dphidr4z; + double dpsidr1x,dpsidr1y,dpsidr1z,dpsidr2x,dpsidr2y,dpsidr2z; + double dpsidr3x,dpsidr3y,dpsidr3z,dpsidr4x,dpsidr4y,dpsidr4z; + + // Definition of cross-term dihedrals + + // phi dihedral + // |--------------------| + // a1-----a2-----a3-----a4-----a5 cross-term atoms + // C N CA C N cross-term atom types + // |--------------------| + // psi dihedral + + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; + + ecmap = 0.0; + int eflag = eflag_caller; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = 0; + + for (n = 0; n < ncrosstermlist; n++) { + i1 = crosstermlist[n][0]; + i2 = crosstermlist[n][1]; + i3 = crosstermlist[n][2]; + i4 = crosstermlist[n][3]; + i5 = crosstermlist[n][4]; + + type = crosstermlist[n][5]; + if (type == 0) continue; + + // calculate bond vectors for both dihedrals + + // phi + // vb21 = r2 - r1 + + vb21x = x[i2][0] - x[i1][0]; + vb21y = x[i2][1] - x[i1][1]; + vb21z = x[i2][2] - x[i1][2]; + vb12x = -1.0*vb21x; + vb12y = -1.0*vb21y; + vb12z = -1.0*vb21z; + vb32x = x[i3][0] - x[i2][0]; + vb32y = x[i3][1] - x[i2][1]; + vb32z = x[i3][2] - x[i2][2]; + vb23x = -1.0*vb32x; + vb23y = -1.0*vb32y; + vb23z = -1.0*vb32z; + + vb34x = x[i3][0] - x[i4][0]; + vb34y = x[i3][1] - x[i4][1]; + vb34z = x[i3][2] - x[i4][2]; + + // psi + // bond vectors same as for phi: vb32 + + vb43x = -1.0*vb34x; + vb43y = -1.0*vb34y; + vb43z = -1.0*vb34z; + + vb45x = x[i4][0] - x[i5][0]; + vb45y = x[i4][1] - x[i5][1]; + vb45z = x[i4][2] - x[i5][2]; + vb54x = -1.0*vb45x; + vb54y = -1.0*vb45y; + vb54z = -1.0*vb45z; + + // calculate normal vectors for planes that define the dihedral angles + + a1x = vb12y*vb23z - vb12z*vb23y; + a1y = vb12z*vb23x - vb12x*vb23z; + a1z = vb12x*vb23y - vb12y*vb23x; + + b1x = vb43y*vb23z - vb43z*vb23y; + b1y = vb43z*vb23x - vb43x*vb23z; + b1z = vb43x*vb23y - vb43y*vb23x; + + a2x = vb23y*vb34z - vb23z*vb34y; + a2y = vb23z*vb34x - vb23x*vb34z; + a2z = vb23x*vb34y - vb23y*vb34x; + + b2x = vb45y*vb43z - vb45z*vb43y; + b2y = vb45z*vb43x - vb45x*vb43z; + b2z = vb45x*vb43y - vb45y*vb43x; + + // calculate terms used later in calculations + + r32 = sqrt(vb32x*vb32x + vb32y*vb32y + vb32z*vb32z); + a1sq = a1x*a1x + a1y*a1y + a1z*a1z; + b1sq = b1x*b1x + b1y*b1y + b1z*b1z; + + r43 = sqrt(vb43x*vb43x + vb43y*vb43y + vb43z*vb43z); + a2sq = a2x*a2x + a2y*a2y + a2z*a2z; + b2sq = b2x*b2x + b2y*b2y + b2z*b2z; + //if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001) + // printf("a1sq b1sq a2sq b2sq: %f %f %f %f \n",a1sq,b1sq,a2sq,b2sq); + if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001) continue; + dpr21r32 = vb21x*vb32x + vb21y*vb32y + vb21z*vb32z; + dpr34r32 = vb34x*vb32x + vb34y*vb32y + vb34z*vb32z; + dpr32r43 = vb32x*vb43x + vb32y*vb43y + vb32z*vb43z; + dpr45r43 = vb45x*vb43x + vb45y*vb43y + vb45z*vb43z; + + // calculate the backbone dihedral angles as VMD and GROMACS + + phi = dihedral_angle_atan2(vb21x,vb21y,vb21z,a1x,a1y,a1z,b1x,b1y,b1z,r32); + psi = dihedral_angle_atan2(vb32x,vb32y,vb32z,a2x,a2y,a2z,b2x,b2y,b2z,r43); + + if (phi == 180.0) phi= -180.0; + if (psi == 180.0) psi= -180.0; + + phi1 = phi; + if (phi1 < 0.0) phi1 += 360.0; + psi1 = psi; + if (psi1 < 0.0) psi1 += 360.0; + + // find the neighbor grid point index + + li1 = int(((phi1+CMAPXMIN2)/CMAPDX)+((CMAPDIM*1.0)/2.0)); + li2 = int(((psi1+CMAPXMIN2)/CMAPDX)+((CMAPDIM*1.0)/2.0)); + + li3 = int((phi-CMAPXMIN2)/CMAPDX); + li4 = int((psi-CMAPXMIN2)/CMAPDX); + mli3 = li3 % CMAPDIM; + mli4 = li4 % CMAPDIM; + mli31 = (li3+1) % CMAPDIM; + mli41 = (li4+1) %CMAPDIM; + mli1 = li1 % CMAPDIM; + mli2 = li2 % CMAPDIM; + mli11 = (li1+1) % CMAPDIM; + mli21 = (li2+1) %CMAPDIM; + t1 = type-1; + if (t1 < 0 || t1 > 5) error->all(FLERR,"Invalid CMAP crossterm_type"); + + // determine the values and derivatives for the grid square points + + gs[0] = cmapgrid[t1][mli3][mli4]; + gs[1] = cmapgrid[t1][mli31][mli4]; + gs[2] = cmapgrid[t1][mli31][mli41]; + gs[3] = cmapgrid[t1][mli3][mli41]; + d1gs[0] = d1cmapgrid[t1][mli1][mli2]; + d1gs[1] = d1cmapgrid[t1][mli11][mli2]; + d1gs[2] = d1cmapgrid[t1][mli11][mli21]; + d1gs[3] = d1cmapgrid[t1][mli1][mli21]; + + d2gs[0] = d2cmapgrid[t1][mli1][mli2]; + d2gs[1] = d2cmapgrid[t1][mli11][mli2]; + d2gs[2] = d2cmapgrid[t1][mli11][mli21]; + d2gs[3] = d2cmapgrid[t1][mli1][mli21]; + + d12gs[0] = d12cmapgrid[t1][mli1][mli2]; + d12gs[1] = d12cmapgrid[t1][mli11][mli2]; + d12gs[2] = d12cmapgrid[t1][mli11][mli21]; + d12gs[3] = d12cmapgrid[t1][mli1][mli21]; + + // calculate the cmap energy and the gradient (dE/dphi,dE/dpsi) + + bc_interpol(phi,psi,li3,li4,gs,d1gs,d2gs,d12gs); + + // sum up cmap energy contributions + + engfraction = 0.2 * E; + if (i1 < nlocal) ecmap += engfraction; + if (i2 < nlocal) ecmap += engfraction; + if (i3 < nlocal) ecmap += engfraction; + if (i4 < nlocal) ecmap += engfraction; + if (i5 < nlocal) ecmap += engfraction; + + // calculate the derivatives dphi/dr_i + + dphidr1x = 1.0*r32/a1sq*a1x; + dphidr1y = 1.0*r32/a1sq*a1y; + dphidr1z = 1.0*r32/a1sq*a1z; + + dphidr2x = -1.0*r32/a1sq*a1x - dpr21r32/a1sq/r32*a1x + + dpr34r32/b1sq/r32*b1x; + dphidr2y = -1.0*r32/a1sq*a1y - dpr21r32/a1sq/r32*a1y + + dpr34r32/b1sq/r32*b1y; + dphidr2z = -1.0*r32/a1sq*a1z - dpr21r32/a1sq/r32*a1z + + dpr34r32/b1sq/r32*b1z; + + dphidr3x = dpr34r32/b1sq/r32*b1x - dpr21r32/a1sq/r32*a1x - r32/b1sq*b1x; + dphidr3y = dpr34r32/b1sq/r32*b1y - dpr21r32/a1sq/r32*a1y - r32/b1sq*b1y; + dphidr3z = dpr34r32/b1sq/r32*b1z - dpr21r32/a1sq/r32*a1z - r32/b1sq*b1z; + + dphidr4x = r32/b1sq*b1x; + dphidr4y = r32/b1sq*b1y; + dphidr4z = r32/b1sq*b1z; + + // calculate the derivatives dpsi/dr_i + + dpsidr1x = 1.0*r43/a2sq*a2x; + dpsidr1y = 1.0*r43/a2sq*a2y; + dpsidr1z = 1.0*r43/a2sq*a2z; + + dpsidr2x = r43/a2sq*a2x + dpr32r43/a2sq/r43*a2x - dpr45r43/b2sq/r43*b2x; + dpsidr2y = r43/a2sq*a2y + dpr32r43/a2sq/r43*a2y - dpr45r43/b2sq/r43*b2y; + dpsidr2z = r43/a2sq*a2z + dpr32r43/a2sq/r43*a2z - dpr45r43/b2sq/r43*b2z; + + dpsidr3x = dpr45r43/b2sq/r43*b2x - dpr32r43/a2sq/r43*a2x - r43/b2sq*b2x; + dpsidr3y = dpr45r43/b2sq/r43*b2y - dpr32r43/a2sq/r43*a2y - r43/b2sq*b2y; + dpsidr3z = dpr45r43/b2sq/r43*b2z - dpr32r43/a2sq/r43*a2z - r43/b2sq*b2z; + + dpsidr4x = r43/b2sq*b2x; + dpsidr4y = r43/b2sq*b2y; + dpsidr4z = r43/b2sq*b2z; + + // calculate forces on cross-term atoms: F = -(dE/dPhi)*(dPhi/dr) + + f1[0] = dEdPhi*dphidr1x; + f1[1] = dEdPhi*dphidr1y; + f1[2] = dEdPhi*dphidr1z; + f2[0] = dEdPhi*dphidr2x + dEdPsi*dpsidr1x; + f2[1] = dEdPhi*dphidr2y + dEdPsi*dpsidr1y; + f2[2] = dEdPhi*dphidr2z + dEdPsi*dpsidr1z; + f3[0] = -dEdPhi*dphidr3x - dEdPsi*dpsidr2x; + f3[1] = -dEdPhi*dphidr3y - dEdPsi*dpsidr2y; + f3[2] = -dEdPhi*dphidr3z - dEdPsi*dpsidr2z; + f4[0] = -dEdPhi*dphidr4x - dEdPsi*dpsidr3x; + f4[1] = -dEdPhi*dphidr4y - dEdPsi*dpsidr3y; + f4[2] = -dEdPhi*dphidr4z - dEdPsi*dpsidr3z; + f5[0] = -dEdPsi*dpsidr4x; + f5[1] = -dEdPsi*dpsidr4y; + f5[2] = -dEdPsi*dpsidr4z; + + // apply force to each of the 5 atoms + + if (i1 < nlocal) { + f[i1][0] += f1[0]; + f[i1][1] += f1[1]; + f[i1][2] += f1[2]; + } + if (i2 < nlocal) { + f[i2][0] += f2[0]; + f[i2][1] += f2[1]; + f[i2][2] += f2[2]; + } + if (i3 < nlocal) { + f[i3][0] += f3[0]; + f[i3][1] += f3[1]; + f[i3][2] += f3[2]; + } + if (i4 < nlocal) { + f[i4][0] += f4[0]; + f[i4][1] += f4[1]; + f[i4][2] += f4[2]; + } + if (i5 < nlocal) { + f[i5][0] += f5[0]; + f[i5][1] += f5[1]; + f[i5][2] += f5[2]; + } + + // tally energy and/or virial + + if (evflag) { + nlist = 0; + if (i1 < nlocal) list[nlist++] = i1; + if (i2 < nlocal) list[nlist++] = i2; + if (i3 < nlocal) list[nlist++] = i3; + if (i4 < nlocal) list[nlist++] = i4; + if (i5 < nlocal) list[nlist++] = i5; + vcmap[0] = (vb12x*f1[0])+(vb32x*f3[0])+((vb43x+vb32x)*f4[0])+ + ((vb54x+vb43x+vb32x)*f5[0]); + vcmap[1] = (vb12y*f1[1])+(vb32y*f3[1])+((vb43y+vb32y)*f4[1])+ + ((vb54y+vb43y+vb32y)*f5[1]); + vcmap[2] = (vb12z*f1[2])+(vb32z*f3[2])+((vb43z+vb32z)*f4[2])+ + ((vb54z+vb43z+vb32z)*f5[2]); + vcmap[3] = (vb12x*f1[1])+(vb32x*f3[1])+((vb43x+vb32x)*f4[1])+ + ((vb54x+vb43x+vb32x)*f5[1]); + vcmap[4] = (vb12x*f1[2])+(vb32x*f3[2])+((vb43x+vb32x)*f4[2])+ + ((vb54x+vb43x+vb32x)*f5[2]); + vcmap[5] = (vb12y*f1[2])+(vb32y*f3[2])+((vb43y+vb32y)*f4[2])+ + ((vb54y+vb43y+vb32y)*f5[2]); + ev_tally(nlist,list,5.0,E,vcmap); + //ev_tally(5,list,nlocal,newton_bond,E,vcmap); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCMAP::post_force_respa(int vflag, int ilevel, int iloop) +{ + if (ilevel == nlevels_respa-1) post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixCMAP::min_post_force(int vflag) +{ + post_force(vflag); +} + +/* ---------------------------------------------------------------------- + energy of CMAP term +------------------------------------------------------------------------- */ + +double FixCMAP::compute_scalar() +{ + double all; + MPI_Allreduce(&ecmap,&all,1,MPI_DOUBLE,MPI_SUM,world); + return all; +} + +// ---------------------------------------------------------------------- +// ---------------------------------------------------------------------- +// methods to read CMAP potential file, perform interpolation +// ---------------------------------------------------------------------- +// ---------------------------------------------------------------------- + +void FixCMAP::read_grid_map(char *cmapfile) +{ + char line[MAXLINE]; + char *chunk; + int i1, i2, i3, i4, i5, i6, j1, j2, j3, j4, j5, j6, counter; + + FILE *fp = fopen(cmapfile,"r"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open fix cmap file %s",cmapfile); + error->one(FLERR,str); + } + + for (int ix1 = 0; ix1 < 6; ix1++) + for (int ix2 = 0; ix2 < CMAPDIM; ix2++) + for (int ix3 = 0; ix3 < CMAPDIM; ix3++) + cmapgrid[ix1][ix2][ix3] = 0.0; + + counter = 0; + i1 = i2 = i3 = i4 = i5 = i6 = 0; + j1 = j2 = j3 = j4 = j5 = j6 = 0; + + while (fgets(line,MAXLINE,fp) != NULL) { + if (line == "" || line[0] == '#') { ;; } + + // read in the cmap grid point values + // NOTE: The order to read the 6 grid maps is HARD-CODED, thus errors + // will occur if content of the file "cmap.data" is altered + // + // Reading order of the maps: + // 1. Alanine map + // 2. Alanine before proline map + // 3. Proline map + // 4. Two adjacent prolines map + // 5. Glycine map + // 6. Glycine before proline map + + else { + chunk = strtok(line, " \r\n"); + while (chunk != NULL) { + + // alanine map + + if (counter < CMAPDIM*CMAPDIM) { + cmapgrid[0][i1][j1] = atof(chunk); + chunk = strtok(NULL, " \r\n"); + j1++; + if (j1 == CMAPDIM) { + j1 = 0; + i1++; + } + counter++; + } + + // alanine-proline map + + else if (counter >= CMAPDIM*CMAPDIM && + counter < 2*CMAPDIM*CMAPDIM) { + cmapgrid[1][i2][j2]= atof(chunk); + chunk = strtok(NULL, " \r\n"); + j2++; + if (j2 == CMAPDIM) { + j2 = 0; + i2++; + } + counter++; + } + + // proline map + + else if (counter >= 2*CMAPDIM*CMAPDIM && + counter < 3*CMAPDIM*CMAPDIM) { + cmapgrid[2][i3][j3] = atof(chunk); + chunk = strtok(NULL, " \r\n"); + j3++; + if (j3 == CMAPDIM) { + j3 = 0; + i3++; + } + counter++; + } + + // 2 adjacent prolines map + + else if (counter >= 3*CMAPDIM*CMAPDIM && + counter < 4*CMAPDIM*CMAPDIM) { + cmapgrid[3][i4][j4] = atof(chunk); + chunk = strtok(NULL, " \r\n"); + j4++; + if (j4 == CMAPDIM) { + j4 = 0; + i4++; + } + counter++; + } + + // glycine map + + else if (counter >= 4*CMAPDIM*CMAPDIM && + counter < 5*CMAPDIM*CMAPDIM) { + cmapgrid[4][i5][j5] = atof(chunk); + chunk = strtok(NULL, " \r\n"); + j5++; + if (j5 == CMAPDIM) { + j5 = 0; + i5++; + } + counter++; + } + + // glycine-proline map + + else if (counter >= 5*CMAPDIM*CMAPDIM && + counter < 6*CMAPDIM*CMAPDIM) { + cmapgrid[5][i6][j6] = atof(chunk); + chunk = strtok(NULL, " \r\n"); + j6++; + if (j6 == CMAPDIM) { + j6 = 0; + i6++; + } + counter++; + } + + else break; + } + } + } + + fclose(fp); +} + +/* ---------------------------------------------------------------------- */ + +void FixCMAP::spline(double *y, double *ddy, int n) +{ + // create the 2nd dervatives of a taublated function y_i(x_i) + // at the tabulated points + + int i, j; + double p, *u; + + memory->create(u,n-1,"cmap:u"); + + ddy[0] = u[0] = 0.0; + + for (i = 1; i <= n-2; i++) { + p = 1.0/(ddy[i-1]+4.0); + ddy[i] = -p; + u[i] = ((((6.0*y[i+1])-(12.0*y[i])+(6.0*y[i-1]))/(CMAPDX*CMAPDX))-u[i-1])*p; + } + + ddy[n-1] = 0.0; + + for (j = n-2; j >= 0; j--) + ddy[j] = ddy[j]*ddy[j+1] + u[j]; + + memory->destroy(u); +} + +/* ---------------------------------------------------------------------- */ + +void FixCMAP::spl_interpolate(double x, double *y, double *ddy, double &yo, + double &dyo) +{ + // perform a 1D cubic spline interpolation + + int ix; + double a,b,a1,b1,a2,b2; + + ix = int((x-CMAPXMIN)/CMAPDX-(1./2.)); + + a = (CMAPXMIN+(ix*1.0)*CMAPDX-x)/CMAPDX; + b = (x-CMAPXMIN-(((ix-1)*1.0)*CMAPDX))/CMAPDX; + + a1 = a*a*a-a; + b1 = b*b*b-b; + + a2 = 3.0*a*a-1.0; + b2 = 3.0*b*b-1.0; + yo = a*y[ix]+b*y[ix+1]+(a1*ddy[ix]+b1*ddy[ix+1])*(CMAPDX*CMAPDX)/6.0; + dyo = (y[ix+1]-y[ix])/CMAPDX-a2/6.0*CMAPDX*ddy[ix]+b2/6.0*CMAPDX*ddy[ix+1]; +} + +/* ---------------------------------------------------------------------- */ + +void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo, + double **d12yo) +{ + // precompute the gradient and cross-derivatives of the map grid points. + // use the bicubic spline to calculate the derivatives + + int i, j, k, ii, jj, xm, p; + double phi, psi, y, d1y, d2y, d12y, tyyk,tdyk; + double *tmp_y, *tmp_dy, *tmp_ddy, **tmap, **tddmap; + int ix; + double a,b,a1,b1,a2,b2; + + xm = CMAPDIM/2; + p = CMAPDIM; + + y = 0.; + d1y = 0.; + d2y = 0.; + d12y = 0.; + + memory->create(tmp_y,CMAPDIM*2,"cmap:tmp_y"); + memory->create(tmp_dy,CMAPDIM*2,"cmap:tmp_dy"); + memory->create(tmp_ddy,CMAPDIM*2,"cmap:tmp_ddy"); + memory->create(tmap,CMAPDIM*2,CMAPDIM*2,"cmap:tmap"); + memory->create(tddmap,CMAPDIM*2,CMAPDIM*2,"cmap:tddmap"); + + // periodically expand the original map + // use the expanded map for bicubic spline interpolation, + // which is used to obtain the derivatives + // actual interpolation is done with bicubic interpolation + + for (i = 0; i < CMAPDIM*2; i++) { + ii = ((i+CMAPDIM-xm)%CMAPDIM); + for (j = 0; j < CMAPDIM*2; j++) { + jj = ((j+CMAPDIM-xm)%CMAPDIM); + tmap[i][j] = map[ii][jj]; + } + } + + for (i = 0; i < CMAPDIM*2; i++) + spline(tmap[i], tddmap[i], CMAPDIM*2); + + for (i = xm; i < CMAPDIM+xm; i++) { + phi = (i-xm)*CMAPDX-180.0; + for (j = xm; j < CMAPDIM+xm; j++) { + psi = (j-xm)*CMAPDX-180.0; + ix = int((psi-CMAPXMIN)/CMAPDX); + a = (CMAPXMIN+((ix+1)*1.0)*CMAPDX-psi)/CMAPDX; + b = (psi-CMAPXMIN-((ix)*1.0)*CMAPDX)/CMAPDX; + a1 = a*a*a-a; + b1 = b*b*b-b; + a2 = 3.0*a*a-1.0; + b2 = 3.0*b*b-1.0; + for (k = 0; k < CMAPDIM*2; k++) { + tyyk = tmp_y[k]; + tdyk = tmp_dy[k]; + tyyk = a*tmap[k][ix]+b*tmap[k][ix+1]+ + (a1*tddmap[k][ix]+b1*tddmap[k][ix+1])*(CMAPDX*CMAPDX)/6.0; + tdyk = (tmap[k][ix+1]-tmap[k][ix])/CMAPDX- + (a2/6.0*CMAPDX*tddmap[k][ix])+(b2/6.0*CMAPDX*tddmap[k][ix+1]); + tmp_y[k] = tyyk; + tmp_dy[k] = tdyk; + } + + spline(tmp_y,tmp_ddy,CMAPDIM+xm+xm); + ix = int((phi-CMAPXMIN)/CMAPDX); + a = (CMAPXMIN+((ix+1)*1.0)*CMAPDX-phi)/CMAPDX; + b = (phi-CMAPXMIN-(ix*1.0)*CMAPDX)/CMAPDX; + a1 = a*a*a-a; + b1 = b*b*b-b; + a2 = 3.0*a*a-1.0; + b2 = 3.0*b*b-1.0; + y = a*tmp_y[ix]+b*tmp_y[ix+1]+ + (a1*tmp_ddy[ix]+b1*tmp_ddy[ix+1])*(CMAPDX*CMAPDX)/6.0; + d1y = (tmp_y[ix+1]-tmp_y[ix])/CMAPDX- + a2/6.0*CMAPDX*tmp_ddy[ix]+b2/6.0*CMAPDX*tmp_ddy[ix+1]; + spline(tmp_dy,tmp_ddy,CMAPDIM+xm+xm); + ix = int((phi-CMAPXMIN)/CMAPDX); + a = (CMAPXMIN+((ix+1)*1.0)*CMAPDX-phi)/CMAPDX; + b = (phi-CMAPXMIN-(ix*1.0)*CMAPDX)/CMAPDX; + a1 = a*a*a-a; + b1 = b*b*b-b; + a2 = 3.0*a*a-1.0; + b2 = 3.0*b*b-1.0; + d2y = a*tmp_dy[ix]+b*tmp_dy[ix+1]+ + (a1*tmp_ddy[ix]+b1*tmp_ddy[ix+1])*(CMAPDX*CMAPDX)/6.0; + d12y = (tmp_dy[ix+1]-tmp_dy[ix])/CMAPDX- + a2/6.0*CMAPDX*tmp_ddy[ix]+b2/6.0*CMAPDX*tmp_ddy[ix+1]; + d1yo[i%p][j%p] = d1y; + d2yo[i%p][j%p] = d2y; + d12yo[i%p][j%p] = d12y; + } + } + + memory->destroy(tmp_y); + memory->destroy(tmp_dy); + memory->destroy(tmp_ddy); + memory->destroy(tmap); + memory->destroy(tddmap); +} + +/* ---------------------------------------------------------------------- */ + +double FixCMAP::dihedral_angle_atan2(double fx, double fy, double fz, + double ax, double ay, double az, + double bx, double by, double bz, + double absg) +{ + // calculate the dihedral angle + + double angle, arg1, arg2; + + arg1 = absg*(fx*bx+fy*by+fz*bz); + arg2 = ax*bx+ay*by+az*bz; + + if (arg1 == 0 && arg2 == 0) + error->all(FLERR,"CMAP: atan2 function cannot take 2 zero arguments"); + else { + angle = atan2(arg1,arg2); + angle = angle*180.0/MY_PI; + } + + return angle; +} + +/* ---------------------------------------------------------------------- */ + +void FixCMAP::bc_coeff(double *gs, double *d1gs, double *d2gs, double *d12gs) +{ + // calculate the bicubic interpolation coefficients c_ij + + static int wt[16][16] = + { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, + -3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, 0, 0, 0, 0, + 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, + 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, + 0, 0, 0, 0,-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, + 0, 0, 0, 0, 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, + -3, 3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0, + 9,-9, 9,-9, 6, 3,-3,-6, 6,-6,-3, 3, 4, 2, 1, 2, + -6, 6,-6, 6,-4,-2, 2, 4,-3, 3, 3,-3,-2,-1,-1,-2, + 2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0, + -6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1, + 4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1 + }; + + int i, j, k, l, in; + double xx, x[16]; + + for (i = 0; i < 4; i++) { + x[i] = gs[i]; + x[i+4] = d1gs[i]*CMAPDX; + x[i+8] = d2gs[i]*CMAPDX; + x[i+12] = d12gs[i]*CMAPDX*CMAPDX; + } + + in = 0; + for (i = 0; i < 4; i++) { + for (j = 0; j < 4; j++) { + xx = 0.0; + for (k = 0; k < 16; k++) xx += wt[in][k]*x[k]; + in++; + cij[i][j] = xx; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixCMAP::bc_interpol(double x1, double x2, int low1, int low2, double *gs, + double *d1gs, double *d2gs, double *d12gs) +{ + // for a given point of interest and its corresponding grid square values, + // gradients and cross-derivatives + // calculate the interpolated value of the point of interest (POI) + + int i, p=12; + double t, u, fac, gs1l, gs2l, gs1u, gs2u; + + // set the interpolation coefficients + + bc_coeff(gs,d1gs,d2gs,d12gs); + + gs1l = g_axis[low1]; + gs2l = g_axis[low2]; + + t = (x1-gs1l)/CMAPDX; + u = (x2-gs2l)/CMAPDX; + + E = dEdPhi = dEdPsi = 0.0; + + for (i = 3; i >= 0; i--) { + E = t*E + ((cij[i][3]*u+cij[i][2])*u+cij[i][1])*u+cij[i][0]; + dEdPhi = u*dEdPhi + (3.0*cij[3][i]*t+2.0*cij[2][i])*t+cij[1][i]; + dEdPsi = t*dEdPsi + (3.0*cij[i][3]*u+2.0*cij[i][2])*u+cij[i][1]; + } + + dEdPhi *= (180.0/MY_PI/CMAPDX); + dEdPsi *= (180.0/MY_PI/CMAPDX); +} + +// ---------------------------------------------------------------------- +// ---------------------------------------------------------------------- +// methods to read and write data file +// ---------------------------------------------------------------------- +// ---------------------------------------------------------------------- + +void FixCMAP::read_data_header(char *line) +{ + if (strstr(line,"crossterms")) { + sscanf(line,BIGINT_FORMAT,&ncmap); + } else error->all(FLERR,"Invalid read data header line for fix cmap"); + + // didn't set in constructor b/c this fix could be defined + // before newton command + + newton_bond = force->newton_bond; +} + +/* ---------------------------------------------------------------------- + unpack N lines in buf from section of data file labeled by keyword + id_offset is applied to atomID fields if multiple data files are read + store CMAP interactions as if newton_bond = OFF, even if actually ON +------------------------------------------------------------------------- */ + +void FixCMAP::read_data_section(char *keyword, int n, char *buf, + tagint id_offset) +{ + int m,tmp,itype; + tagint atom1,atom2,atom3,atom4,atom5; + char *next; + + next = strchr(buf,'\n'); + *next = '\0'; + int nwords = atom->count_words(buf); + *next = '\n'; + + if (nwords != 7) { + char str[128]; + sprintf(str,"Incorrect %s format in data file",keyword); + error->all(FLERR,str); + } + + // loop over lines of CMAP crossterms + // tokenize the line into values + // add crossterm to one of my atoms, depending on newton_bond + + for (int i = 0; i < n; i++) { + next = strchr(buf,'\n'); + *next = '\0'; + sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT + " " TAGINT_FORMAT " " TAGINT_FORMAT, + &tmp,&itype,&atom1,&atom2,&atom3,&atom4,&atom5); + + atom1 += id_offset; + atom2 += id_offset; + atom3 += id_offset; + atom4 += id_offset; + atom5 += id_offset; + + if ((m = atom->map(atom1)) >= 0) { + if (num_crossterm[m] == CMAPMAX) + error->one(FLERR,"Too many CMAP crossterms for one atom"); + crossterm_type[m][num_crossterm[m]] = itype; + crossterm_atom1[m][num_crossterm[m]] = atom1; + crossterm_atom2[m][num_crossterm[m]] = atom2; + crossterm_atom3[m][num_crossterm[m]] = atom3; + crossterm_atom4[m][num_crossterm[m]] = atom4; + crossterm_atom5[m][num_crossterm[m]] = atom5; + num_crossterm[m]++; + } + + if ((m = atom->map(atom2)) >= 0) { + if (num_crossterm[m] == CMAPMAX) + error->one(FLERR,"Too many CMAP crossterms for one atom"); + crossterm_type[m][num_crossterm[m]] = itype; + crossterm_atom1[m][num_crossterm[m]] = atom1; + crossterm_atom2[m][num_crossterm[m]] = atom2; + crossterm_atom3[m][num_crossterm[m]] = atom3; + crossterm_atom4[m][num_crossterm[m]] = atom4; + crossterm_atom5[m][num_crossterm[m]] = atom5; + num_crossterm[m]++; + } + + if ((m = atom->map(atom3)) >= 0) { + if (num_crossterm[m] == CMAPMAX) + error->one(FLERR,"Too many CMAP crossterms for one atom"); + crossterm_type[m][num_crossterm[m]] = itype; + crossterm_atom1[m][num_crossterm[m]] = atom1; + crossterm_atom2[m][num_crossterm[m]] = atom2; + crossterm_atom3[m][num_crossterm[m]] = atom3; + crossterm_atom4[m][num_crossterm[m]] = atom4; + crossterm_atom5[m][num_crossterm[m]] = atom5; + num_crossterm[m]++; + } + + if ((m = atom->map(atom4)) >= 0) { + if (num_crossterm[m] == CMAPMAX) + error->one(FLERR,"Too many CMAP crossterms for one atom"); + crossterm_type[m][num_crossterm[m]] = itype; + crossterm_atom1[m][num_crossterm[m]] = atom1; + crossterm_atom2[m][num_crossterm[m]] = atom2; + crossterm_atom3[m][num_crossterm[m]] = atom3; + crossterm_atom4[m][num_crossterm[m]] = atom4; + crossterm_atom5[m][num_crossterm[m]] = atom5; + num_crossterm[m]++; + } + + if ((m = atom->map(atom5)) >= 0) { + if (num_crossterm[m] == CMAPMAX) + error->one(FLERR,"Too many CMAP crossterms for one atom"); + crossterm_type[m][num_crossterm[m]] = itype; + crossterm_atom1[m][num_crossterm[m]] = atom1; + crossterm_atom2[m][num_crossterm[m]] = atom2; + crossterm_atom3[m][num_crossterm[m]] = atom3; + crossterm_atom4[m][num_crossterm[m]] = atom4; + crossterm_atom5[m][num_crossterm[m]] = atom5; + num_crossterm[m]++; + } + + buf = next + 1; + } +} + +/* ---------------------------------------------------------------------- */ + +bigint FixCMAP::read_data_skip_lines(char *keyword) +{ + return ncmap; +} + +/* ---------------------------------------------------------------------- + write Mth header line to file + only called by proc 0 +------------------------------------------------------------------------- */ + +void FixCMAP::write_data_header(FILE *fp, int mth) +{ + fprintf(fp,BIGINT_FORMAT " cmap crossterms\n",ncmap); +} + +/* ---------------------------------------------------------------------- + return size I own for Mth data section + # of data sections = 1 for this fix + nx = # of crossterms owned by my local atoms + if newton_bond off, atom only owns crossterm if it is atom3 + ny = columns = type + 5 atom IDs +------------------------------------------------------------------------- */ + +void FixCMAP::write_data_section_size(int mth, int &nx, int &ny) +{ + int i,m; + + tagint *tag = atom->tag; + int nlocal = atom->nlocal; + + nx = 0; + for (i = 0; i < nlocal; i++) + for (m = 0; m < num_crossterm[i]; m++) + if (crossterm_atom3[i][m] == tag[i]) nx++; + + ny = 6; +} + +/* ---------------------------------------------------------------------- + pack values for Mth data section into 2d buf + buf allocated by caller as owned crossterms by 6 +------------------------------------------------------------------------- */ + +void FixCMAP::write_data_section_pack(int mth, double **buf) +{ + int i,m; + + // 1st column = CMAP type + // 2nd-6th columns = 5 atom IDs + + tagint *tag = atom->tag; + int nlocal = atom->nlocal; + + int n = 0; + for (i = 0; i < nlocal; i++) { + for (m = 0; m < num_crossterm[i]; m++) { + if (crossterm_atom3[i][m] != tag[i]) continue; + buf[n][0] = ubuf(crossterm_type[i][m]).d; + buf[n][1] = ubuf(crossterm_atom1[i][m]).d; + buf[n][2] = ubuf(crossterm_atom2[i][m]).d; + buf[n][3] = ubuf(crossterm_atom3[i][m]).d; + buf[n][4] = ubuf(crossterm_atom4[i][m]).d; + buf[n][5] = ubuf(crossterm_atom5[i][m]).d; + n++; + } + } +} + +/* ---------------------------------------------------------------------- + write section keyword for Mth data section to file + use Molecules or Charges if that is only field, else use fix ID + only called by proc 0 +------------------------------------------------------------------------- */ + +void FixCMAP::write_data_section_keyword(int mth, FILE *fp) +{ + fprintf(fp,"\nCMAP\n\n"); +} + +/* ---------------------------------------------------------------------- + write N lines from buf to file + convert buf fields to int or double depending on styles + index can be used to prepend global numbering + only called by proc 0 +------------------------------------------------------------------------- */ + +void FixCMAP::write_data_section(int mth, FILE *fp, + int n, double **buf, int index) +{ + for (int i = 0; i < n; i++) + fprintf(fp,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT + " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT "\n", + index+i,(int) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i, + (tagint) ubuf(buf[i][2]).i,(tagint) ubuf(buf[i][3]).i, + (tagint) ubuf(buf[i][4]).i,(tagint) ubuf(buf[i][5]).i); +} + +// ---------------------------------------------------------------------- +// ---------------------------------------------------------------------- +// methods for restart and communication +// ---------------------------------------------------------------------- +// ---------------------------------------------------------------------- + +/* ---------------------------------------------------------------------- + pack entire state of Fix into one write +------------------------------------------------------------------------- */ + +void FixCMAP::write_restart(FILE *fp) +{ + if (comm->me == 0) { + int size = sizeof(bigint); + fwrite(&size,sizeof(int),1,fp); + fwrite(&ncmap,sizeof(bigint),1,fp); + } +} + +/* ---------------------------------------------------------------------- + use state info from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixCMAP::restart(char *buf) +{ + ncmap = *((bigint *) buf); +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based arrays for restart file +------------------------------------------------------------------------- */ + +int FixCMAP::pack_restart(int i, double *buf) +{ + int n = 1; + for (int m = 0; m < num_crossterm[i]; m++) { + buf[n++] = ubuf(MAX(crossterm_type[i][m],-crossterm_type[i][m])).d; + buf[n++] = ubuf(crossterm_atom1[i][m]).d; + buf[n++] = ubuf(crossterm_atom2[i][m]).d; + buf[n++] = ubuf(crossterm_atom3[i][m]).d; + buf[n++] = ubuf(crossterm_atom4[i][m]).d; + buf[n++] = ubuf(crossterm_atom5[i][m]).d; + } + buf[0] = n; + + return n; +} + +/* ---------------------------------------------------------------------- + unpack values from atom->extra array to restart the fix +------------------------------------------------------------------------- */ + +void FixCMAP::unpack_restart(int nlocal, int nth) +{ + double **extra = atom->extra; + + // skip to Nth set of extra values + + int n = 0; + for (int i = 0; i < nth; i++) n += static_cast (extra[nlocal][n]); + + int count = static_cast (extra[nlocal][n++]); + num_crossterm[nlocal] = (count-1)/6; + + for (int m = 0; m < num_crossterm[nlocal]; m++) { + crossterm_type[nlocal][m] = (int) ubuf(extra[nlocal][n++]).i; + crossterm_atom1[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i; + crossterm_atom2[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i; + crossterm_atom3[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i; + crossterm_atom4[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i; + crossterm_atom5[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i; + } +} + +/* ---------------------------------------------------------------------- + maxsize of any atom's restart data +------------------------------------------------------------------------- */ + +int FixCMAP::maxsize_restart() +{ + return 1 + CMAPMAX*6; +} + +/* ---------------------------------------------------------------------- + size of atom nlocal's restart data +------------------------------------------------------------------------- */ + +int FixCMAP::size_restart(int nlocal) +{ + return 1 + num_crossterm[nlocal]*6; +} + +/* ---------------------------------------------------------------------- + allocate atom-based array +------------------------------------------------------------------------- */ + +void FixCMAP::grow_arrays(int nmax) +{ + num_crossterm = memory->grow(num_crossterm,nmax,"cmap:num_crossterm"); + crossterm_type = memory->grow(crossterm_type,nmax,CMAPMAX, + "cmap:crossterm_type"); + crossterm_atom1 = memory->grow(crossterm_atom1,nmax,CMAPMAX, + "cmap:crossterm_atom1"); + crossterm_atom2 = memory->grow(crossterm_atom2,nmax,CMAPMAX, + "cmap:crossterm_atom2"); + crossterm_atom3 = memory->grow(crossterm_atom3,nmax,CMAPMAX, + "cmap:crossterm_atom3"); + crossterm_atom4 = memory->grow(crossterm_atom4,nmax,CMAPMAX, + "cmap:crossterm_atom4"); + crossterm_atom5 = memory->grow(crossterm_atom5,nmax,CMAPMAX, + "cmap:crossterm_atom5"); + + // must initialize num_crossterm to 0 for added atoms + // may never be set for some atoms when data file is read + + for (int i = nmax_previous; i < nmax; i++) num_crossterm[i] = 0; + nmax_previous = nmax; +} + +/* ---------------------------------------------------------------------- + copy values within local atom-based array +------------------------------------------------------------------------- */ + +void FixCMAP::copy_arrays(int i, int j, int delflag) +{ + num_crossterm[j] = num_crossterm[i]; + + for (int k = 0; k < num_crossterm[j]; k++){ + crossterm_type[j][k] = crossterm_type[i][k]; + crossterm_atom1[j][k] = crossterm_atom1[i][k]; + crossterm_atom2[j][k] = crossterm_atom2[i][k]; + crossterm_atom3[j][k] = crossterm_atom3[i][k]; + crossterm_atom4[j][k] = crossterm_atom4[i][k]; + crossterm_atom5[j][k] = crossterm_atom5[i][k]; + } +} + +/* ---------------------------------------------------------------------- + initialize one atom's array values, called when atom is created +------------------------------------------------------------------------- */ + +void FixCMAP::set_arrays(int i) +{ + num_crossterm[i] = 0; +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based array for exchange with another proc +------------------------------------------------------------------------- */ + +int FixCMAP::pack_exchange(int i, double *buf) +{ + int n = 0; + buf[n++] = ubuf(num_crossterm[i]).d; + for (int m = 0; m < num_crossterm[i]; m++) { + buf[n++] = ubuf(crossterm_type[i][m]).d; + buf[n++] = ubuf(crossterm_atom1[i][m]).d; + buf[n++] = ubuf(crossterm_atom2[i][m]).d; + buf[n++] = ubuf(crossterm_atom3[i][m]).d; + buf[n++] = ubuf(crossterm_atom4[i][m]).d; + buf[n++] = ubuf(crossterm_atom5[i][m]).d; + } + return n; +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based array from exchange with another proc +------------------------------------------------------------------------- */ + +int FixCMAP::unpack_exchange(int nlocal, double *buf) +{ + int n = 0; + num_crossterm[nlocal] = (int) ubuf(buf[n++]).i; + for (int m = 0; m < num_crossterm[nlocal]; m++) { + crossterm_type[nlocal][m] = (int) ubuf(buf[n++]).i; + crossterm_atom1[nlocal][m] = (tagint) ubuf(buf[n++]).i; + crossterm_atom2[nlocal][m] = (tagint) ubuf(buf[n++]).i; + crossterm_atom3[nlocal][m] = (tagint) ubuf(buf[n++]).i; + crossterm_atom4[nlocal][m] = (tagint) ubuf(buf[n++]).i; + crossterm_atom5[nlocal][m] = (tagint) ubuf(buf[n++]).i; + } + return n; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double FixCMAP::memory_usage() +{ + int nmax = atom->nmax; + double bytes = nmax * sizeof(int); // num_crossterm + bytes += nmax*CMAPMAX * sizeof(int); // crossterm_type + bytes += 5*nmax*CMAPMAX * sizeof(int); // crossterm_atom 12345 + bytes += maxcrossterm*6 * sizeof(int); // crosstermlist + return bytes; +} diff --git a/src/MOLECULE/fix_cmap.h b/src/MOLECULE/fix_cmap.h new file mode 100644 index 0000000000..213d881e23 --- /dev/null +++ b/src/MOLECULE/fix_cmap.h @@ -0,0 +1,129 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(cmap,FixCMAP) + +#else + +#ifndef LMP_FIX_CMAP_H +#define LMP_FIX_CMAP_H + +#include "fix.h" +namespace LAMMPS_NS { + +class FixCMAP : public Fix { + public: + FixCMAP(class LAMMPS *, int, char **); + ~FixCMAP(); + int setmask(); + void init(); + void setup(int); + void setup_pre_neighbor(); + void min_setup(int); + void pre_neighbor(); + void pre_reverse(int, int); + void post_force(int); + void post_force_respa(int, int, int); + void min_post_force(int); + double compute_scalar(); + + void read_data_header(char *); + void read_data_section(char *, int, char *, tagint); + bigint read_data_skip_lines(char *); + void write_data_header(FILE *, int); + void write_data_section_size(int, int &, int &); + void write_data_section_pack(int, double **); + void write_data_section_keyword(int, FILE *); + void write_data_section(int, FILE *, int, double **, int); + + void write_restart(FILE *); + void restart(char *); + int pack_restart(int, double *); + void unpack_restart(int, int); + int size_restart(int); + int maxsize_restart(); + + void grow_arrays(int); + void copy_arrays(int, int, int); + void set_arrays(int); + int pack_exchange(int, double *); + int unpack_exchange(int, double *); + + double memory_usage(); + + private: + int nprocs,me; + int newton_bond,eflag_caller; + int ctype,nlevels_respa; + int ncrosstermtypes,crossterm_per_atom,maxcrossterm; + int ncrosstermlist; + bigint ncmap; + + int **crosstermlist; + + int nmax_previous; + int *num_crossterm; + int **crossterm_type; + tagint **crossterm_atom1,**crossterm_atom2,**crossterm_atom3; + tagint **crossterm_atom4,**crossterm_atom5; + + double E,dEdPhi,dEdPsi; + double ecmap; + double fcmap[4],cij[4][4]; + double *g_axis; + + // CMAP grid points obtained from external file + + double ***cmapgrid; + + // partial derivatives and cross-derivatives of the grid data + + double ***d1cmapgrid,***d2cmapgrid,***d12cmapgrid; + + // read map grid data + + void read_grid_map(char *); + + // read in CMAP cross terms from LAMMPS data file + + void read_cmap_data(int, char *); + + // pre-compute the partial and cross-derivatives of map grid points + + void set_map_derivatives(double **, double **, double **, double **); + + // cubic spline interpolation functions for derivatives of map grid points + + void spline(double *, double *, int); + void spl_interpolate(double, double *, double *, double &, double &); + + // calculate dihedral angles + + double dihedral_angle_atan2(double, double, double, double, double, double, + double, double, double, double); + + // calculate bicubic interpolation coefficient matrix c_ij + + void bc_coeff(double *, double *, double *, double *); + + // perform bicubic interpolation at point of interest + + void bc_interpol(double, double, int, int, double *, double *, double *, + double *); +}; +} + +#endif +#endif diff --git a/src/version.h b/src/version.h index 3f90b4fc5b..35c0f59e69 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "21 Sep 2016" +#define LAMMPS_VERSION "26 Sep 2016" diff --git a/tools/ch2lmp/README b/tools/ch2lmp/README index eb994ccb6f..fa031bd757 100755 --- a/tools/ch2lmp/README +++ b/tools/ch2lmp/README @@ -11,35 +11,40 @@ In this directory, you should find: 3) An "example" folder containing an example of how to use these tools. 4) An "other" folder containing other potentially useful tools. -In addition, you will need to provide the following input for +In addition, you will need to provide the following input for charmm2lammps.pl: + 1) a CHARMM parameter file (par_.prm) 2) a CHARMM topology file (top_.rtf) 3) a CHARMM coordinates file or pdb file (.crd or .pdb) + (NOTE: a .pdb file is required if the -cmap option is used) 4) a CHARMM psf file (.psf) To use the charmm2lammps.pl script, type: "perl charmm2lammps.pl -[-option[=#]] " where is the name -of the CHARMM FF you are using, (i.e. all22_prot), and is -the common name of your *.crd and *.psf files. The options for the -script are listed below. If the option requires a parameter, the -syntax must be [-option[=#]], (i.e. -border=5). + [-option=value ...]" where is the +name of the CHARMM FF you are using, (i.e. all22_prot), and +is the common name of your *.crd and *.psf files. Possible options +are listed next. If the option requires a parameter, the syntax +should be -option=value (e.g. -border=5). --help "display available options", --charmm "add charmm types to LAMMPS data file", --water "add TIP3P water [default: 1 g/cc]", --ions "add (counter)ions using Na+ and Cl- [default: 0 mol/l]", --center "recenter atoms", --quiet "do not print info", --pdb_ctrl "output project_ctrl.pdb", --l "set x-, y-, and z-dimensions simultaneously", --lx "x-dimension of simulation box", --ly "y-dimension of simulation box", --lz "z-dimension of simulation box", --border "add border to all sides of simulation box [default: 0 A]", --ax "rotation around x-axis", --ay "rotation around y-axis", --az "rotation around z-axis" +-help display available options +-charmm add charmm types to LAMMPS data file +-water add TIP3P water [default: 1 g/cc] +-ions add (counter)ions using Na+ and Cl- [default: 0 mol/l] +-center recenter atoms +-quiet do not print info +-pdb_ctrl output project_ctrl.pdb +-l set x-, y-, and z-dimensions simultaneously +-lx x-dimension of simulation box +-ly y-dimension of simulation box +-lz z-dimension of simulation box +-border add border to all sides of simulation box [default: 0 A] +-ax rotation around x-axis +-ay rotation around y-axis +-az rotation around z-axis +-cd correction for dihedral for carbohydrate systems +-cmap add CMAP section to data file and fix cmap command lines in + input script" (NOTE: requires use of *.pdb file) In the "example" folder, you will find example files that were created by following the steps below. These steps describe how to take a @@ -100,7 +105,22 @@ molecule. The -pdb_ctrl option produces the 1ac7_ctrl.pdb file that can be visualized in a standard visualization package such as VMD. The -charmm option put comments into the LAMMPS data file (everything after the # sign is a comment) for user convenience in tracking atom -types etc. according to CHARMM nomenclature. +types etc. according to CHARMM nomenclature. + +The example molecule provided above (i.e., 1ac7) is a DNA fragment. +If instead, a peptide longer than 2 amino acid residues or a protein +is to be modeled, the '-cmap' option should be used. This will add a +section at the end of the data file with the heading of 'CMAP' that +will contain cmap crossterm corrections for the phi-psi dihedrals for +the amino acid residues. You will then need to also copy the +appropriate file for the cmap crossterms into your directory and be +sure that you are using the appropriate cmap crossterms that go with +the respective version of the charmm force field that is being used +(e.g, cmap22.data or cmap36.data). This is necessary to account for +the fact that the CHARMM group has provided updated cmap correction +terms for use with the c36 and more recent version of the charmm +protein force field. Copies of cmap22.data and cmap36.data are +provided in the tools/ch2lmp directory. The default timestep in the LAMMPS *.in file is set to 0.5 fs, which can typically be increased to 2 fs after equilibration if the bonds @@ -109,9 +129,9 @@ delay on neigh_modify can probably increased to 5 or so to improve speed. The -ions option allows the user to neutralize the simulation cell -with Na+ or Cl- counterions if the molecule has a net -charge. Additional salt can be added by increasing the default -concentration (i.e. -ions=0.5). +with Na+ or Cl- counterions if the molecule has a net charge +Additional salt can be added by increasing the default concentration +(e.g. -ions=0.5). ** In the "other" file folder, you will find: @@ -126,5 +146,3 @@ concentration (i.e. -ions=0.5). 3) A 3rd party perl script called "crd2pdb.pl" 4) A 3rd party fortran code called "pdb_to_crd.f" - - diff --git a/tools/ch2lmp/charmm2lammps.pl b/tools/ch2lmp/charmm2lammps.pl index 8e43a06cab..0e39f589da 100644 --- a/tools/ch2lmp/charmm2lammps.pl +++ b/tools/ch2lmp/charmm2lammps.pl @@ -28,22 +28,43 @@ # 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. +# 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 # # 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. +# ----------------------------------------------------------------------------- + # Initialization sub Test { - my $name = shift(@_); + my $name = shift(@_); # "@_" = argument passed to the subroutine printf("Error: file %s not found\n", $name) if (!scalar(stat($name))); return !scalar(stat($name)); @@ -54,11 +75,13 @@ { my $k = 0; my @dir = ("x", "y", "z"); - my @options = ("-help", "-nohints", "-water", "-ions", "-center", + + # Modified by Xiaohu Hu, May 2009. Options "-cmap" and "-cdihedral" added + my @options = ("-help", "-charmm", "-water", "-ions", "-center", "-quiet", "-pdb_ctrl", "-l", "-lx", "-ly", "-lz", - "-border", "-ax", "-ay", "-az"); - my @remarks = ("display this message", - "do not print type and style hints in data file", + "-border", "-ax", "-ay", "-az", "-cd", "-cmap"); + my @remarks = ("display this message", + "add charmm types to LAMMPS data file", "add TIP3P water [default: 1 g/cc]", "add (counter)ions using Na+ and Cl- [default: 0 mol/l]", "recenter atoms", @@ -71,14 +94,20 @@ "add border to all sides of simulation box [default: 0 A]", "rotation around x-axis", "rotation around y-axis", - "rotation around z-axis" + "rotation around z-axis", + "use the 6-ring dihedral correction", + "generate the CMAP section" ); my $notes; $program = "charmm2lammps"; - $version = "1.8.3"; - $year = "2016"; - $add = 1; +# $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; $water_dens = 0; $ions = 0; $info = 1; @@ -112,7 +141,7 @@ foreach (@ARGV) { - if (substr($_, 0, 1) eq "-") + if (substr($_, 0, 1) eq "-") # if the first letter of the aguement is "-" = an option { my $k = 0; my @tmp = split("="); @@ -123,23 +152,24 @@ 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 - print("Warning: ignoring unknown command line flag: $tmp[0]\n"); + $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 } else { @@ -150,7 +180,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\n\n", + 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", $program, $version, $year); printf("Usage:\n %s.pl [-option[=#] ..] forcefield project\n\n",$program); printf("Options:\n"); @@ -161,7 +191,7 @@ 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("%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") @@ -300,7 +330,7 @@ sub PSFConnectivity { - my $n = PSFGoto(bonds); + my $n = PSFGoto(bonds); # $n = the total number of bonds return if (scalar(@nconnect)); printf("Info: creating connectivity\n") if ($info); @@ -312,7 +342,6 @@ } } - sub PSFDihedrals # hack to accomodate { # LAMMPS' way of calc $idihedral = 0; # LJ 1-4 interactions @@ -343,7 +372,6 @@ return $ndihedral; } - sub CreatePSFIndex # make an index of id { # locations my @psf_ids = ("!NATOM","!NBOND:","!NTHETA:","!NPHI:","!NIMPHI:"); @@ -363,26 +391,24 @@ } } - - sub PSFGoto # goto $ident in + sub PSFGoto # goto $ident in and return the total number of $ident { CreatePSFIndex() if (!scalar(%PSFIndex)); my $id = shift(@_); - my @n = split(" ", $PSFIndex{$id}); + my @n = split(" ", $PSFIndex{$id}); # = 1 if the $ident is found @PSFBuffer = (); # return PSFDihedrals() if ($id eq "dihedrals"); if (!scalar(@n)) { printf("Warning: PSF index for $id not found\n"); - seek(PSF, 0, SEEK_END); + seek(PSF, 0, SEEK_END); # set file-handle position to the EOF return -1; } seek(PSF, $n[0], SEEK_SET); return $n[1]; } - sub PSFGet { if ($dihedral_flag) @@ -400,7 +426,7 @@ } - sub PSFWrite + sub PSFWrite { my $items = shift(@_); my $n = $items; @@ -516,14 +542,11 @@ sub Markers { - my %markers = ( - NONBONDED => '0', - BONDS => '1', - ANGLES => '2', - DIHEDRALS => '3', - IMPROPERS => '4', - IMPROPER => '4' - ); + my %markers; + my $n = 0; + + foreach ("NONBONDED", "BONDS", "ANGLES", "DIHEDRALS", "IMPROPER") { + $markers{$_} = $n++; } return %markers; } @@ -613,7 +636,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; @@ -688,7 +711,7 @@ sub SetScreeningFactor # set screening factor { - my $id = shift(@_); + my $id = shift(@_); # @_ has the form: ($variable, float no.) my $value = shift(@_); my $new = ""; @@ -702,7 +725,6 @@ $parms[$id] = $new; } - sub CorrectDihedralParameters { my $n = PSFGoto(dihedrals); @@ -713,15 +735,20 @@ my $first; my $last; - for (my $i=0; $i<$n; ++$i) + for (my $i=0; $i<$n; ++$i) # loop over all dihedrals { - 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 @@ -1010,7 +1037,7 @@ 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 "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); @@ -1114,7 +1141,7 @@ { 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 %11.8g %11.8g %11.8g%s\n", + 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 ". @@ -1164,7 +1191,7 @@ CRDGoto(atoms); $net_charge = 0; - printf(LAMMPS "Atoms # full\n\n") if ($n>0); + printf(LAMMPS "Atoms\n\n") if ($n>0); for (my $i=0; $i<$n; ++$i) { my @crd = $pdb ? NextPDB2CRD() : split(" ", ); @@ -1172,7 +1199,7 @@ 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 %11.7g %11.7g %11.7g%s\n", ++$k, + 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 ". @@ -1198,7 +1225,6 @@ { 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; @@ -1212,7 +1238,7 @@ @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\n\n", $header); 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); @@ -1428,15 +1449,29 @@ 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 "pair_style lj/charmm/coul/long 8 10\n"); + 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"); - printf(LAMMPS "read_data $project.data\n\n"); # read data + + 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"); + }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 + } + if ($coefficients ne "") # corrected coeffs { foreach (split(":", $coefficients)) @@ -1446,16 +1481,15 @@ 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"); @@ -1463,6 +1497,972 @@ 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 # + # # + # --------------------------------------------------------------------------- # + +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 +{ + 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 = ; + 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 = ; + 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> <#2>' \n"; + print "where '<#1> is the next atom number, 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"; + +# End of the CharmmCmap subroutine +} # main @@ -1470,4 +2470,6 @@ WriteData(); WriteLAMMPSInput(); printf("Info: conversion complete\n\n") if ($info); + DihedralCorrect6Ring() if ($cdihedral); + CharmmCmap() if ($cmap);