From 8d63e582804e754546d69f2e944c71e534115cb7 Mon Sep 17 00:00:00 2001
From: James Graham <J.A.Graham@soton.ac.uk>
Date: Wed, 17 Aug 2016 15:50:17 +0100
Subject: [PATCH] Clean up forcefield output code

---
 data/martini_v2.2.itp  | 26 +---------------
 data/w.itp             | 30 ++++++++++++++++++
 pycgtool/forcefield.py | 71 ++++++++++++++++++------------------------
 pycgtool/pycgtool.py   |  2 +-
 4 files changed, 62 insertions(+), 67 deletions(-)

diff --git a/data/martini_v2.2.itp b/data/martini_v2.2.itp
index 4f0ba77..fa16da9 100755
--- a/data/martini_v2.2.itp
+++ b/data/martini_v2.2.itp
@@ -907,28 +907,4 @@ BP4 72.0 0.000 A 0.0 0.0
   Qa 	SQ0 	1 	0.19402E-00 	0.20914E-02 ; almost attractive
   SQa 	Q0 	1 	0.19402E-00 	0.20914E-02 ; almost attractive
   SQa 	SQ0 	1 	0.85338E-01 	0.53946E-03 ; 75almost attractive, s=0.43
-  Q0 	SQ0 	1 	0.15091E-00 	0.16267E-02 ; intermediate 
-
-
-  
-;;;;;; WATER (representing 4 H2O molecules)
-
-[ moleculetype ]
-; molname  	nrexcl
-  W 	    	1
-
-[ atoms ]
-;id 	type 	resnr 	residu 	atom 	cgnr 	charge
- 1 	P4 	1 	W 	W 	1 	0 
-
-;;;;;; ANTIFREEZE (prevents freezing of water)
-
-[ moleculetype ]
-; molname 	 nrexcl
-  WF	   	 1
-                                                                                
-[ atoms ]
-;id 	type 	resnr 	residu 	atom 	cgnr 	charge
- 1 	BP4 	1 	WF	WF	1 	0
-
-
+  Q0 	SQ0 	1 	0.15091E-00 	0.16267E-02 ; intermediate
diff --git a/data/w.itp b/data/w.itp
index dbf1bf8..f11253e 100644
--- a/data/w.itp
+++ b/data/w.itp
@@ -1,3 +1,33 @@
+; MARTINI FORCEFIELD V2.2
+;
+; SJ MARRINK (last modified: 04-12-2012 by DdJ)
+;
+; NOTE 1: Bead definitinions in this file have not been changed with respect to
+;         martini_v2.1.itp. This file is purely here for clarity sake. Differences
+;         between V2.1 and V2.2 are created by the martinize script and can be found
+;         in martini_v2.2_aminoacids.itp.
+;
+; please cite:
+;
+; D.H. de Jong, G. Singh, W.F.D. Bennet, C. Arnarez, T.A. Wassenaar, L.V. Schafer,
+; X. Periole, D.P. Tieleman, S.J. Marrink.
+; Improved Parameters for the Martini Coarse-Grained Protein Force Field
+; J. Chem. Theory Comput., DOI: 10.1021/ct300646g
+
+; L. Monticelli, S. Kandasamy, X. Periole, R. Larson, D.P. Tieleman, S.J. Marrink.
+; The MARTINI coarse grained force field: extension to proteins.
+; J. Chem. Th. Comp., 4:819-834, 2008.
+;
+; S.J. Marrink, H.J. Risselada, S. Yefimov, D.P. Tieleman, A.H. de Vries.
+; The MARTINI forcefield: coarse grained model for biomolecular simulations.
+; JPC-B, 111:7812-7824, 2007.
+;
+; and (if using lipid topologies):
+;
+; S.J. Marrink, A.H. de Vries, A.E. Mark.
+; Coarse grained model for semi-quantitative lipid simulations.
+; JPC-B, 108:750-760, 2004.
+
 ;;;;;; WATER (representing 4 H2O molecules)
 
 [ moleculetype ]
diff --git a/pycgtool/forcefield.py b/pycgtool/forcefield.py
index 5f83b49..2b4190d 100644
--- a/pycgtool/forcefield.py
+++ b/pycgtool/forcefield.py
@@ -13,36 +13,38 @@ class ForceField:
     """
     Class used to output a GROMACS .ff forcefield
     """
+    # TODO output r2b file for ends of chain
     def __init__(self, name):
         """
-        Open a forcefield directory at name.  If it does not exist it is created.
+        Open a named forcefield directory.  If it does not exist it is created.
 
-        :param name: Directory name to open/create
+        :param name: Forcefield name to open/create
         """
-        os.makedirs(name, exist_ok=True)
+        self.dirname = "ff{0}.ff".format(name)
+        os.makedirs(self.dirname, exist_ok=True)
 
-        self.dirname = name
         with open(os.path.join(self.dirname, "forcefield.itp"), "w") as itp:
-            print("#define _FF_PYCGTOOL", file=itp)
+            print("#define _FF_PYCGTOOL_{0}".format(name), file=itp)
             print('#include "martini_v2.2.itp"', file=itp)
 
+        dist_dat_dir = os.path.join(dir_up(os.path.realpath(__file__), 2), "data")
         # Copy main MARTINI itp
-        shutil.copyfile(os.path.join(dir_up(os.path.realpath(__file__), 2), "data", "martini_v2.2.itp"),
+        shutil.copyfile(os.path.join(dist_dat_dir, "martini_v2.2.itp"),
                         os.path.join(self.dirname, "martini_v2.2.itp"))
         # Copy water models
-        shutil.copyfile(os.path.join(dir_up(os.path.realpath(__file__), 2), "data", "watermodels.dat"),
+        shutil.copyfile(os.path.join(dist_dat_dir, "watermodels.dat"),
                         os.path.join(self.dirname, "watermodels.dat"))
-        shutil.copyfile(os.path.join(dir_up(os.path.realpath(__file__), 2), "data", "w.itp"),
+        shutil.copyfile(os.path.join(dist_dat_dir, "w.itp"),
                         os.path.join(self.dirname, "w.itp"))
 
         # Create atomtypes.atp required for correct masses with pdb2gmx
-        with CFG(os.path.join(dir_up(os.path.realpath(__file__), 2), "data", "martini_v2.2.itp"), allow_duplicate=True) as cfg,\
+        with CFG(os.path.join(dist_dat_dir, "martini_v2.2.itp"), allow_duplicate=True) as cfg,\
                 open(os.path.join(self.dirname, "atomtypes.atp"), 'w') as atomtypes:
             for toks in cfg["atomtypes"]:
                 print(" ".join(toks), file=atomtypes)
 
         with open(os.path.join(self.dirname, "forcefield.doc"), "w") as doc:
-            print("PyCGTOOL produced MARTINI force field", file=doc)
+            print("PyCGTOOL produced MARTINI force field - {0}".format(name), file=doc)
 
     def write_rtp(self, name, mapping, bonds):
         """
@@ -54,50 +56,37 @@ class ForceField:
         :param mapping: AA->CG mapping from which to collect molecules
         :param bonds: BondSet from which to collect bonds
         """
-        # TODO print everything to file
+        def write_bond_angle_dih(bonds, section_header, file, multiplicity=None):
+            if bonds:
+                print("  [ {0:s} ]".format(section_header), file=file)
+            for bond in bonds:
+                line = "    " + " ".join(["{0:>4s}".format(atom) for atom in bond.atoms])
+                line += " {0:12.5f} {1:12.5f}".format(bond.eqm, bond.fconst)
+                if multiplicity is not None:
+                    line += " {0:4d}".format(multiplicity)
+                print(line, file=file)
+
         with open(os.path.join(self.dirname, name), "w") as rtp:
             print("[ bondedtypes ]", file=rtp)
             print(("{:4d}" * 8).format(1, 1, 1, 1, 1, 1, 0, 0), file=rtp)
 
             for mol in mapping:
-                try:
-                    bonds[mol]
-                except KeyError:
-                    # Skip molecules without bonds
+                # Skip molecules without bonds
+                if mol not in bonds:
                     continue
 
                 print("[ {0} ]".format(mol), file=rtp)
 
-                print(" [ atoms ]", file=rtp)
+                print("  [ atoms ]", file=rtp)
                 for bead in mapping[mol]:
-                    #        name  type  charge  c-group
-                    print("  {:4s} {:4s} {:3.6f} {:4d}".format(
+                    #          name  type  charge  chg-group
+                    print("    {:>4s} {:>4s} {:3.6f} {:4d}".format(
                         bead.name, bead.type, bead.charge, 0
                     ), file=rtp)
 
                 bond_tmp = bonds.get_bond_lengths(mol, with_constr=True)
-                if bond_tmp:
-                    print(" [ bonds ]", file=rtp)
-                for bond in bond_tmp:
-                    print("  {:4s} {:4s} {:8.3f} {:8.3f}".format(
-                        bond.atoms[0], bond.atoms[1],
-                        bond.eqm, bond.fconst
-                    ), file=rtp)
-
+                write_bond_angle_dih(bond_tmp, "bonds", rtp)
                 bond_tmp = bonds.get_bond_angles(mol)
-                if bond_tmp:
-                    print(" [ angles ]", file=rtp)
-                for bond in bond_tmp:
-                    print("  {:4s} {:4s} {:4s} {:8.3f} {:8.3f}".format(
-                        bond.atoms[0], bond.atoms[1], bond.atoms[2],
-                        bond.eqm, bond.fconst
-                    ), file=rtp)
-
+                write_bond_angle_dih(bond_tmp, "angles", rtp)
                 bond_tmp = bonds.get_bond_dihedrals(mol)
-                if bond_tmp:
-                    print(" [ dihedrals ]", file=rtp)
-                for bond in bond_tmp:
-                    print("  {:4s} {:4s} {:4s} {:4s} {:8.3f} {:8.3f} {:4d}".format(
-                        bond.atoms[0], bond.atoms[1], bond.atoms[2], bond.atoms[3],
-                        bond.eqm, bond.fconst, 1
-                    ), file=rtp)
+                write_bond_angle_dih(bond_tmp, "dihedrals", rtp, multiplicity=1)
diff --git a/pycgtool/pycgtool.py b/pycgtool/pycgtool.py
index f0f4473..0654167 100755
--- a/pycgtool/pycgtool.py
+++ b/pycgtool/pycgtool.py
@@ -67,7 +67,7 @@ def main(args, config):
             print("If your CG molecule should contain charges the itp will need to be edited")
             if config.output_forcefield:
                 logger.info("Creating GROMACS forcefield directory")
-                ff = ForceField("ff" + config.output_name + ".ff")
+                ff = ForceField(config.output_name)
                 ff.write_rtp(config.output_name + ".rtp", mapping, bonds)
                 logger.info("GROMACS forcefield directory created")
             else:
-- 
GitLab