diff --git a/1-lj-liquid.py b/1-lj-liquid.py new file mode 100644 index 0000000000000000000000000000000000000000..7ad22f0b67df42a2eb8533ba66e090c3e729ef80 --- /dev/null +++ b/1-lj-liquid.py @@ -0,0 +1,107 @@ + +# 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 +)