Coverage for examples/molecules_polydisperse.py: 100%
104 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4Created on Thu Mar 13 13:33:58 2025
6@author: nbailey
7"""
9import numpy as np
10import gamdpy as gp
12gp.select_gpu()
14# Simulation params
15rho, temperature = 0.85, 1.5
16N_nominal = 2000
17chain_lengths = [5, 10]
18composition = [1, 2] # non-negative integers. Should be relatively prime (no common factors)
20# size_base is the number of atoms in the smallest repeating unit
21size_base = 0
22for cl, frac in zip(chain_lengths, composition):
23 size_base += cl*frac
25num_base_units = N_nominal // size_base
26num_mols_each_type = []
27N_mol = 0
28N = 0
30for cl, frac in zip(chain_lengths, composition):
31 num_mols_each_type.append(frac*num_base_units)
32 N_mol += num_mols_each_type[-1]
33 N += num_mols_each_type[-1] * cl
35#N_mol = N_nominal//chain_length
36#N = N_mol * chain_length
39print(f"N={N}; N_mol={N_mol};num_base_units={num_base_units}")
40print("num_mols_each_type")
41print(num_mols_each_type)
43filename = 'Data/chains_poly'
44num_timeblocks_equilibration = 64
45num_timeblocks_production = 16
46steps_per_timeblock = 1 * 1024
48molecule_dicts = []
51for index, cl in enumerate(chain_lengths):
52 pos_this_mol = []
53 types_this_mol = []
54 masses_this_mol = []
55 for i in range(cl):
56 pos_this_mol.append( [ i*1.0, (i%2)*.1, 0. ] ) # x, y, z for this particle
57 types_this_mol.append( index ) # all particles in a molecule have the same type, but this differs from molecule to molecule
58 masses_this_mol.append( 1.0 )
60 # Setup configuration: Single molecule first, then duplicate
61 top_this_mol = gp.Topology([f'MyMolecule{cl}', ])
62 top_this_mol.bonds = gp.bonds_from_positions(pos_this_mol, cut_off=1.1, bond_type=0)
63 top_this_mol.angles = gp.angles_from_bonds(top_this_mol.bonds, angle_type=0)
64 top_this_mol.dihedrals = gp.dihedrals_from_angles(top_this_mol.angles, dihedral_type=0)
65 top_this_mol.molecules[f'MyMolecule{cl}'] = gp.molecules_from_bonds(top_this_mol.bonds)
68 dict_this_mol = {"positions" : pos_this_mol,
69 "particle_types" : types_this_mol,
70 "masses" : masses_this_mol,
71 "topology" : top_this_mol}
72 molecule_dicts.append(dict_this_mol)
74 print(f'Initial Positions for molecule with chain length: {cl}')
75 for position in pos_this_mol:
76 print('\t\t', position)
77 print('Particle types:\t', types_this_mol)
78 print('Bonds: \t', top_this_mol.bonds)
79 print('Angles: \t', top_this_mol.angles)
80 print('Dihedrals: \t', top_this_mol.dihedrals)
81 print()
84 # This call creates the pdf "molecule.pdf" with a drawing of the molecule
85 # Use block=True to visualize the molecule before running the simulation
86 gp.plot_molecule(top_this_mol, pos_this_mol, types_this_mol, filename=f"chain{cl}.pdf", block=False)
89configuration = gp.replicate_molecules(molecule_dicts, num_mols_each_type, safety_distance=3.0)
90configuration.randomize_velocities(temperature=temperature)
92print(f'Number of molecules: {len(configuration.topology.molecules[f"MyMolecule{chain_lengths[0]}"])}, {len(configuration.topology.molecules[f"MyMolecule{chain_lengths[1]}"])}')
93print(f'Number of particles: {configuration.N}\n')
95# Make bond interactions
96bond_potential = gp.harmonic_bond_function
97bond_params = [[0.8, 1000.], ]
98bonds = gp.Bonds(bond_potential, bond_params, configuration.topology.bonds)
100# Make angle interactions
101angle0, k = 2.0, 500.0
102angle_potential = gp.cos_angle_function
103angles = gp.Angles(angle_potential, configuration.topology.angles, parameters=[[k, angle0],])
105# Make dihedral interactions
106rbcoef=[.0, 5.0, .0, .0, .0, .0]
107dihedral_potential = gp.ryckbell_dihedral
108dihedrals = gp.Dihedrals(dihedral_potential, configuration.topology.dihedrals, parameters=[rbcoef, ])
110# Exlusion list
111exclusions = dihedrals.get_exclusions(configuration)
112#exclusions = bonds.get_exclusions(configuration)
114# Make pair potential
115pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon)
116sig = [[1.00, 0.95],
117 [0.95, 0.9]]
118eps = [[1.00, 0.95],
119 [0.95, 0.90]]
120cut = np.array(sig)*2.5
124pair_pot = gp.PairPotential(pair_func, params=[sig, eps, cut], exclusions=exclusions, max_num_nbs=1000)
126# Make integrator
127integrator = gp.integrators.NVT(temperature=temperature, tau=0.1, dt=0.004)
129# Setup runtime actions, i.e. actions performed during simulation of timeblocks
130runtime_actions = [gp.TrajectorySaver(),
131 gp.ScalarSaver(),
132 gp.MomentumReset(100)]
134# Setup simulation
135sim = gp.Simulation(configuration, [pair_pot, bonds, angles, dihedrals], integrator, runtime_actions,
136 num_timeblocks=num_timeblocks_equilibration, steps_per_timeblock=steps_per_timeblock,
137 storage='memory')
139print('\nCompression and equilibration: ')
140dump_filename = 'Data/dump_compress.lammps'
141with open(dump_filename, 'w') as f:
142 print(gp.configuration_to_lammps(sim.configuration, timestep=0), file=f)
145initial_rho = configuration.N / configuration.get_volume()
146for block in sim.run_timeblocks():
147 volume = configuration.get_volume()
148 N = configuration.N
149 current_rho = N/volume
150 print(sim.status(per_particle=True), f'rho= {current_rho:.3}', end='\t')
151 print(f'P= {(N*temperature + np.sum(configuration["W"]))/volume:.3}') # pV = NkT + W
152 with open(dump_filename, 'a') as f:
153 print(gp.configuration_to_lammps(sim.configuration, timestep=sim.steps_per_block*(block+1)), file=f)
155 # Scale configuration to get closer to final density, rho
156 if block<sim.num_blocks/2:
157 desired_rho = (block+1)/(sim.num_blocks/2)*(rho - initial_rho) + initial_rho
158 if desired_rho > 1.2*current_rho:
159 desired_rho = 1.2*current_rho
160 configuration.atomic_scale(density=desired_rho)
161 configuration.copy_to_device() # Since we altered configuration, we need to copy it back to device
162print(sim.summary())
163print(configuration)
165sim = gp.Simulation(configuration, [pair_pot, bonds, angles, dihedrals], integrator, runtime_actions,
166 num_timeblocks=num_timeblocks_production, steps_per_timeblock=steps_per_timeblock,
167 compute_plan=sim.compute_plan, storage=filename+'.h5')
169print('\nProduction: ')
170dump_filename = 'Data/dump.lammps'
171with open(dump_filename, 'w') as f:
172 print(gp.configuration_to_lammps(sim.configuration, timestep=0), file=f)
175for block in sim.run_timeblocks():
176 print(sim.status(per_particle=True))
177 with open(dump_filename, 'a') as f:
178 print(gp.configuration_to_lammps(sim.configuration, timestep=sim.steps_per_block*(block+1)), file=f)
180print(sim.summary())
181print(configuration)