Merge pull request #22 from polysimtools/release-0.2.2

Release 0.2.2

pysimm/cassandra.py

@@ -718,6 +718,13 @@ def __init__(self, sst, **kwargs):
                                      'Assuming it is **Type-1** force field')
 
                     sst.ff_class = '1'
+
+                if not all([pt.name for pt in sst.particle_types]):
+                    self.logger.error('The name of at least one particle type in MC system is not defined. '
+                                      'Will not be able to map particles back after the CASSANDRA simulations. '
+                                      '\nPlease, setup the names for all particle types for your MC system')
+                    exit(1)
+
             # Decorating the system with bonds_fixed flag and angle_fixed flag
             for bt in sst.bond_types:
                 bt.is_fixed = True
@@ -1275,6 +1282,7 @@ class McfWriter(object):
     def __init__(self, syst, file_ref):
         self.syst = syst
         self.file_ref = file_ref
+        self.logger = logging.getLogger('MCF Writer')
 
     def write(self, typing='all'):
         """pysimm.cassandra.McfWriter.write
@@ -1315,26 +1323,32 @@ def __write_atom_info__(self, out):
             count = 0
             line_template = '{l[0]:<6}{l[1]:<7}{l[2]:<5}{l[3]:<8.3f}{l[4]:<10.6f}' \
                             '{l[5]:<6}{l[6]:<11.3f}{l[7]:<9.3f}\n'
+            warn_flag = False
             for item in self.syst.particles:
-                line = [count + 1, '', '', '', 0, 'LJ', 0, 0]
-                if hasattr(item, 'charge'):
+                line = [count + 1, '', '', 0, 0, 'LJ', 0, 0]
+                if item.charge:
                     line[4] = item.charge
-                else:
-                    line[4] = 0
-                if hasattr(item, 'type'):
-                    if hasattr(item.type, 'name'):
+                if item.type:
+                    line[1] = item.type.tag
+                    line[2] = item.type.tag
+                    if item.type.name:
                         line[1] = item.type.name
-                    if hasattr(item.type, 'elem'):
                         line[2] = item.type.elem
-                    if hasattr(item.type, 'mass'):
+                    else:
+                        warn_flag = True
+                    if item.type.mass:
                         line[3] = item.type.mass
-                    if hasattr(item.type, 'epsilon'):
+                    if item.type.epsilon:
                         line[6] = KCALMOL_2_K * item.type.epsilon
+                    if item.type.sigma:
                         line[7] = item.type.sigma
                 else:
                     continue
                 out.write(line_template.format(l=line))
                 count += 1
+            if warn_flag:
+                self.logger.warning('Some particle type names (and/or element names) inside the system are not defined.'
+                                    ' Will use type identifiers instead')
         else:
             self.__write_empty__(out, text_tag)
         out.write('\n')

tests/test_cassandra.py

@@ -1,41 +1,60 @@
 import unittest
 from os import path as osp
 from pysimm import system, cassandra
+import pytest
 
 
 class CassandraTestCase(unittest.TestCase):
 
+    data_path = 'test_data'
+
     def test_gas_unwrap(self):
 
         # Read molecular system from the file and decorate the frame  with fixed property
-        test_sst = system.read_lammps(osp.join('test_data', 'wrapped_co2.lmps'))
+        test_sst = system.read_lammps(osp.join(self.data_path, 'wrapped_co2.lmps'))
         css = cassandra.Cassandra(test_sst)
         for p in css.system.particles:
-            if p.tag > 3500:
-                p.is_fixed = False
+            p.is_fixed = False
 
         # Assert that gas molecules are wrapped
         bnd_lng = []
-        for ind in [3501, 3504, 3507]:
-            for i in [1, 2]:
-                bnd_lng.append(
-                    (css.system.particles[ind].x - css.system.particles[ind + i].x) ** 2 +
-                    (css.system.particles[ind].y - css.system.particles[ind + i].y) ** 2 +
-                    (css.system.particles[ind].z - css.system.particles[ind + i].z) ** 2)
+        for p in test_sst.particles:
+            if p.type.name == 'C':
+                for i in [1, 2]:
+                    bnd_lng.append(
+                        (css.system.particles[p.tag].x - css.system.particles[p.tag + i].x) ** 2 +
+                        (css.system.particles[p.tag].y - css.system.particles[p.tag + i].y) ** 2 +
+                        (css.system.particles[p.tag].z - css.system.particles[p.tag + i].z) ** 2)
         self.assertFalse(all([(1.344 < el) and (el < 1.346) for el in bnd_lng]))
 
         css.unwrap_gas()
 
         # Assert that gas molecules are unwrapped
         bnd_lng = []
-        for ind in [3501, 3504, 3507]:
-            for i in [1, 2]:
-                bnd_lng.append(
-                    (css.system.particles[ind].x - css.system.particles[ind + i].x) ** 2 +
-                    (css.system.particles[ind].y - css.system.particles[ind + i].y) ** 2 +
-                    (css.system.particles[ind].z - css.system.particles[ind + i].z) ** 2)
+        for p in test_sst.particles:
+            if p.type.name == 'C':
+                for i in [1, 2]:
+                    bnd_lng.append(
+                        (css.system.particles[p.tag].x - css.system.particles[p.tag + i].x) ** 2 +
+                        (css.system.particles[p.tag].y - css.system.particles[p.tag + i].y) ** 2 +
+                        (css.system.particles[p.tag].z - css.system.particles[p.tag + i].z) ** 2)
         self.assertTrue(all([(1.344 < el) and (el < 1.346) for el in bnd_lng]))
 
+    def test_nonames_mc(self):
+        sst = system.System()
+        sst.dim = system.Dimension(dx=40, dy=40, dz=40, center=[0, 0, 0])
+        sst.forcefield = 'trappe/amber'
+        css = cassandra.Cassandra(sst)
+        my_gcmc_props = css.read_input(osp.join(self.data_path, 'props.inp'))
+
+        specie = system.read_lammps(osp.join(self.data_path, 'toluene_nonames.lmps'))
+        specie.forcefield = 'trappe/amber'
+
+        with pytest.raises(SystemExit) as err_back:
+            css.add_gcmc(species=specie, is_rigid=True, max_ins=200,
+                         chem_pot=-30.34, out_folder=self.data_path, **my_gcmc_props)
+            assert err_back.type == SystemExit
+
 
 if __name__ == '__main__':
     unittest.main()

tests/test_data/props.inp

@@ -0,0 +1,47 @@
+! Basic template input file for pysimm.Cassandra for GCMC in MC-MD iterations
+!   it lists "static' properties that will be set once for all MC simulations
+!   during MC-MD. It is made for convinience. Alternatively all this properties
+!   might be set from the MC-MD Python script.
+
+# VDW_Style
+lj cut_tail 14.0
+
+# Rcutoff_Low
+0.0
+
+# Mixing_Rule
+lb
+
+# Charge_Style
+coul ewald 14.0 1e-5
+
+# Temperature_Info
+300
+
+# Run_Type
+Equilibration 100
+
+# Simulation_Length_Info
+Units       Steps
+prop_freq   500
+coord_freq  10000
+run         10000
+
+# Average_Info
+1
+
+# Property_Info 1
+Energy_Total
+Density
+Nmols
+Pressure
+
+# CBMC_Info
+kappa_ins 12
+kappa_dih 10
+rcut_cbmc 2.0
+
+# Pair_Energy
+true
+
+END

tests/test_data/toluene_nonames.lmps

@@ -0,0 +1,78 @@
+read using pysimm.system.read_mol
+
+7 atoms
+7 bonds
+8 angles
+0 dihedrals
+0 impropers
+
+3 atom types
+2 bond types
+1 angle types
+
+-13.290500 13.292500 xlo xhi
+-12.121700 12.841700 ylo yhi
+-10.549200 11.028200 zlo zhi
+
+Masses
+
+   1	15.034
+   2	13.018
+   3	12.01
+
+Pair Coeffs
+
+   1	0.19474	3.75
+   2	0.10035	3.695
+   3	0.04173	3.88
+
+Bond Coeffs
+
+   1	5000.0	1.4
+   2	5000.0	1.54
+
+Angle Coeffs
+
+   1	5000.0	120.0
+
+Atoms
+
+   1	1	3	0	-1.21241	-0.331716	0.000460468
+   2	1	2	0	1.21246	-0.331491	-0.000693376
+   3	1	2	0	9.25687e-05	-1.0316	1.03773e-05
+   4	1	2	0	-1.21254	1.06828	0.000398264
+   5	1	2	0	1.21233	1.06851	-0.000389329
+   6	1	2	0	-0.000167649	1.7684	3.28067e-07
+   7	1	1	0	-2.54602	-1.10184	-0.000427083
+
+Velocities
+
+   1	-1.92247e-07	-1.50612e-07	-1.49082e-08
+   2	1.94278e-07	-1.51941e-07	5.54101e-08
+   3	4.35049e-09	2.42428e-08	-1.47441e-08
+   4	-1.76794e-07	-3.05153e-09	7.56444e-09
+   5	1.7681e-07	-3.67578e-09	-2.01054e-08
+   6	-1.57793e-09	1.63953e-07	-5.34978e-10
+   7	2.02747e-07	4.76605e-08	3.77879e-09
+
+Bonds
+
+   1	1	1	3
+   2	1	1	4
+   3	2	1	7
+   4	1	2	3
+   5	1	2	5
+   6	1	4	6
+   7	1	5	6
+
+Angles
+
+   1	1	3	1	4
+   2	1	3	1	7
+   3	1	4	1	7
+   4	1	3	2	5
+   5	1	1	3	2
+   6	1	6	4	1
+   7	1	6	5	2
+   8	1	4	6	5
+