Add old code file

1-lj-liquid.py

+
+# coding: utf-8
+
+# # Basic usage of the `arbdmodel` package
+
+# ### Step 1: Create particle types 
+
+# In[1]:
+
+
+import numpy as np
+from arbdmodel import ParticleType, PointParticle, ArbdModel, ArbdEngine
+
+system_size = 620
+num_particles = 6400
+
+## diffusion of liquid argon: https://www.sciencedirect.com/science/article/pii/0375960171901290
+## units "60 * 3.41 AA * sqrt(119.5 k K / 40 dalton)" "AA**2/ns"
+argon = ParticleType(name="Ar",
+                     mass=39.95,              # dalton
+                     damping_coefficient=5000, # per ns
+                     custom_property="my_value", 
+                     epsilon=0.177,  # kcal_mol
+                     radius=3.345/2) # Angstroms
+
+print(argon)
+
+
+# ### Step 2: Create a system
+
+# In[2]:
+
+
+## Generate Nx3 array of random cartesian coordinates
+positions = system_size*(np.random.random([num_particles,3])-0.5)
+
+## Create a list of point particles located at those positions
+particles = [PointParticle(type_=argon, position=positions[i])
+                        for i in range(num_particles)]
+
+model = ArbdModel(particles, 
+                  dimensions=[system_size for i in range(3)], # Ångstroms
+                  timestep = 2e-6                             # ns; can be specified below with engine instead
+                 )
+print(model)
+print(model.children[:2])
+
+
+# ### Step 3: Describe the interactions between the particles
+# 
+
+# In[3]:
+
+
+from arbdmodel.interactions import LennardJones
+
+lj = LennardJones()
+model.add_nonbonded_interaction(lj)
+
+
+# ### Step 4: Run the simulation
+
+# In[4]:
+
+engine = ArbdEngine(output_period = 1e2, num_steps = 1e5)
+engine.simulate(model, output_name="1-lj", directory='sims/1-argon',
+                num_steps=1e3, gpu=1
+)                               # simulation parameters here override those in constructor
+
+
+### Step 5: Customize the interactions
+from arbdmodel.interactions import NonbondedPotential
+
+## Clear nonbonded interactions
+model.nonbonded_interactions = []
+
+class BuckinghamPotential(NonbondedPotential):
+    def potential(self, r, types):
+        ## https://royalsocietypublishing.org/doi/10.1098/rspa.1938.0173
+        ## optionally could pull information from typeA, typeB
+        typeA, typeB = types
+        return 143932.65 * ( 1.69 * np.exp( -r/0.273 ) - 102e-4 / r**6 )
+
+pot = BuckinghamPotential()
+model.add_nonbonded_interaction( pot )
+
+engine.simulate(model, output_name="2-buck", directory='sims/1-argon',
+                num_steps = 1e3, gpu=1
+)
+
+### Step 6: Add bonds
+from arbdmodel.interactions import HarmonicBond
+bond = HarmonicBond(k=1,           # kcal/mol AA**2
+                    r0 = 3         # AA
+)
+
+## Go through every other index
+for i in range(len(0,model.children,2):
+    j = i+1
+    model.add_bond( model.children[i],
+                    model.children[j],
+                    bond )
+print(len(model.bonds))
+
+engine.simulate(model, output_name="3-bonds", directory='sims/1-argon',
+                num_steps = 1e3, gpu=1
+)