Merge pull request #51 from ICAMS/fix_mass

Fix mass

Author: srmnitc

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.2.3
+current_version = 1.2.4
 commit = True
 tag = True
 

calphy/__init__.py

@@ -4,7 +4,7 @@
 from calphy.alchemy import Alchemy
 from calphy.routines import MeltingTemp
 
-__version__ = "1.2.3"
+__version__ = "1.2.4"
 
 def addtest(a,b):
     return a+b

calphy/alchemy.py

@@ -81,7 +81,7 @@ def run_averaging(self):
         lmp = ph.create_structure(lmp, self.calc)
 
         #set up potential
-        lmp = ph.set_potential(lmp, self.calc)
+        lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
 
         #add some computes
         lmp.command("variable         mvol equal vol")
@@ -149,7 +149,7 @@ def run_integration(self, iteration=1):
 
         #set up hybrid potential
         #here we only need to set one potential
-        lmp = ph.set_potential(lmp, self.calc)
+        lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
         #lmp = ph.set_double_hybrid_potential(lmp, self.options, self.calc._pressureair_style, self.calc._pressureair_coeff)
 
         #remap the box to get the correct pressure

calphy/helpers.py

@@ -99,7 +99,19 @@ def create_structure(lmp, calc, species=None):
     return lmp
 
 
-def set_potential(lmp, options):
+def set_mass(lmp, options, ghost_elements=0):
+    count = 1
+    for i in range(options.n_elements):
+        lmp.mass(i+1, options.mass[i])
+        count += 1
+
+    for i in range(ghost_elements):
+        lmp.mass(count+i, 1.00)
+
+    return lmp
+
+
+def set_potential(lmp, options, ghost_elements=0):
     """
     Set the interatomic potential
 
@@ -116,8 +128,8 @@ def set_potential(lmp, options):
     lmp.pair_style(options.pair_style[0])
     lmp.pair_coeff(options.pair_coeff[0])
 
-    for i in range(options.n_elements):
-        lmp.mass(i+1, options.mass[i])
+    lmp = set_mass(lmp, options, ghost_elements=ghost_elements)
+
     return lmp
 
 
@@ -158,7 +170,7 @@ def get_structures(file, species, index=None):
         aseobjs = traj[index].to_ase(species=species)
     return aseobjs
 
-def set_hybrid_potential(lmp, options, eps):
+def set_hybrid_potential(lmp, options, eps, ghost_elements=0):
     pc =  options.pair_coeff[0]
     pcraw = pc.split()
     pcnew = " ".join([*pcraw[:2], *[options.pair_style[0],], *pcraw[2:]])
@@ -167,11 +179,11 @@ def set_hybrid_potential(lmp, options, eps):
     lmp.command("pair_coeff       %s"%pcnew)
     lmp.command("pair_coeff       * * ufm %f 1.5"%eps) 
 
-    for i in range(options.n_elements):
-        lmp.mass(i+1, options.mass[i])
+    lmp = set_mass(lmp, options, ghost_elements=ghost_elements)
+
     return lmp
 
-def set_double_hybrid_potential(lmp, options, pair_style, pair_coeff):
+def set_double_hybrid_potential(lmp, options, pair_style, pair_coeff, ghost_elements=0):
 
     pc1 =  pair_coeff[0]
     pcraw1 = pc1.split()
@@ -191,8 +203,8 @@ def set_double_hybrid_potential(lmp, options, pair_style, pair_coeff):
     lmp.command("pair_coeff       %s"%pcnew1)
     lmp.command("pair_coeff       %s"%pcnew2) 
 
-    for i in range(options.n_elements):
-        lmp.mass(i+1, options.mass[i])
+    lmp = set_mass(lmp, options, ghost_elements=ghost_elements)
+
     return lmp
 
 def remap_box(lmp, x, y, z):

calphy/input.py

@@ -209,7 +209,7 @@ def __init__(self):
 
         self.tolerance = InputTemplate()
         self.tolerance.lattice_constant = 0.0002
-        self.tolerance.spring_constant = 0.01
+        self.tolerance.spring_constant = 0.1
         self.tolerance.solid_fraction = 0.7
         self.tolerance.liquid_fraction = 0.05
         self.tolerance.pressure = 0.5

calphy/liquid.py

@@ -118,7 +118,7 @@ def run_averaging(self):
         lmp = ph.create_structure(lmp, self.calc, species=self.calc.n_elements+self.calc._ghost_element_count)
 
         #set up potential
-        lmp = ph.set_potential(lmp, self.calc)
+        lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
 
         #Melt regime for the liquid
         lmp.velocity("all create", self.calc._temperature_high, np.random.randint(0, 10000))
@@ -183,7 +183,7 @@ def run_integration(self, iteration=1):
 
         #set hybrid ufm and normal potential
         #lmp = ph.set_hybrid_potential(lmp, self.options, self.eps)
-        lmp = ph.set_potential(lmp, self.calc)
+        lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
 
         #remap the box to get the correct pressure
         lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)

calphy/phase.py

@@ -675,7 +675,7 @@ def reversible_scaling(self, iteration=1):
         lmp = ph.read_dump(lmp, conf, species=self.calc.n_elements+self.calc._ghost_element_count)
 
         #set up potential
-        lmp = ph.set_potential(lmp, self.calc)
+        lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
 
         #remap the box to get the correct pressure
         lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)
@@ -755,7 +755,7 @@ def reversible_scaling(self, iteration=1):
         else:
             self.check_if_solidfied(lmp, "traj.temp.dat")
 
-        lmp = ph.set_potential(lmp, self.calc)
+        lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
 
         #reverse scaling
         lmp.command("variable         flambda equal ramp(${li},${lf})")
@@ -850,7 +850,7 @@ def temperature_scaling(self, iteration=1):
         lmp = ph.read_dump(lmp, conf, species=self.calc.n_elements+self.calc._ghost_element_count)
 
         #set up potential
-        lmp = ph.set_potential(lmp, self.calc)
+        lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
 
         #remap the box to get the correct pressure
         lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)
@@ -936,7 +936,7 @@ def pressure_scaling(self, iteration=1):
         lmp = ph.read_dump(lmp, conf, species=self.calc.n_elements+self.calc._ghost_element_count)
 
         #set up potential
-        lmp = ph.set_potential(lmp, self.calc)
+        lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
 
         #remap the box to get the correct pressure
         lmp = ph.remap_box(lmp, self.lx, self.ly, self.lz)

calphy/solid.py

@@ -143,7 +143,7 @@ def run_averaging(self):
 
         #set up potential
         if self.calc.potential_file is None:
-            lmp = ph.set_potential(lmp, self.calc)
+            lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
         else:
             lmp.command("include %s"%self.calc.potential_file)
 
@@ -214,7 +214,7 @@ def run_integration(self, iteration=1):
 
         #set up potential
         if self.calc.potential_file is None:
-            lmp = ph.set_potential(lmp, self.calc)
+            lmp = ph.set_potential(lmp, self.calc, ghost_elements=self.calc._ghost_element_count)
         else:
             lmp.command("include %s"%self.calc.potential_file)
 

setup.py

@@ -60,7 +60,7 @@
     test_suite='tests',
     tests_require=test_requirements,
     url='https://github.com/ICAMS/calphy',
-    version='1.2.3',
+    version='1.2.4',
     zip_safe=False,
     entry_points={
         'console_scripts': [