Coverage for examples/tethered_particles.py: 97%

1""" Example using tethered LJ particles  

2 

3 Particles with label 1 and 2 are tethered with different Hooke springs.  

4 All particles are integrated forward in time with the NVT integrator. 

5 Density of fluid/free particles and tethered particles is the same. 

6 Initial and final configurations are saved in xyz format for easy inspection in vmd  

7 Both ways of defining the anchor points (from particle indices and types) 

8 are shown. 

9""" 

10import os 

11import numpy as np 

12import gamdpy as gp 

13 

14# Setup a default fcc configuration 

15use_list = True 

16 

17nx, ny, nz, rho = 6, 6, 10, 1.0 

18configuration = gp.Configuration(D=3) 

19configuration.make_lattice(gp.unit_cells.FCC, cells=[nx, ny, nz], rho=rho) 

20configuration['m'] = 1.0 

21configuration.randomize_velocities(temperature=2.0) 

22 

23# Fluid/free particles have type '0', tethered particles '1' and '2' 

24indices = []; kspring = []; 

25for n in range(configuration.N): 

26 if -3 < configuration['r'][n][2] < -2: 

27 indices.append(n) 

28 kspring.append(300) 

29 configuration.ptype[n] = 1 

30 elif -1 < configuration['r'][n][1] < 0: 

31 indices.append(n) 

32 kspring.append(500) 

33 configuration.ptype[n] = 2 

34 

35gp.tools.save_configuration(configuration, "initial.xyz") 

36 

37# Tether specifications.  

38tether = gp.Tether() 

39 

40if use_list: 

41 tether.set_anchor_points_from_lists(indices, kspring, configuration) 

42else: 

43 tether.set_anchor_points_from_types(particle_types=[1, 2], spring_constants=[300, 500], configuration=configuration) 

44 

45# Set the pair interactions 

46pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

47sig = [[1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, 1.0, 1.0]] 

48eps = [[1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, 1.0, 1.0]] 

49cut = np.array(sig)*2.5 

50pair_pot = gp.PairPotential(pair_func, params=[sig, eps, cut], max_num_nbs=1000) 

51 

52# Setup integrator: NVT 

53integrator = gp.integrators.NVT(temperature=2.0, tau=0.2, dt=0.005) 

54 

55# Setup runtime actions, i.e. actions performed during simulation of timeblocks 

56runtime_actions = [gp.TrajectorySaver(), 

57 gp.ScalarSaver()] 

58 

59# Setup Simulation. Total number of time steps: num_blocks * steps_per_block 

60sim = gp.Simulation(configuration, [pair_pot, tether], integrator, runtime_actions, 

61 num_timeblocks=16, steps_per_timeblock=1024, 

62 storage='memory') 

63 

64# Run simulation one block at a time 

65for block in sim.run_timeblocks(): 

66 print(sim.status(per_particle=True)) 

67 

68print(sim.summary()) 

69 

70gp.tools.save_configuration(configuration, "final.xyz") 

71