Coverage for examples/poiseuille.py: 100%

1""" Example of a nanoslit pore simulation using tethered LJ particles  

2 

3 The particles are tethered with a Hooke spring. The wall particles interact with a relaxation device.  

4 All particles are integrated forward in time with the NVE integrator 

5 

6 Wall density is set to 1.0 and fluid density to 0.8 - this is achieved by 

7 inclusion of a dummy particle type  

8 

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

10""" 

11 

12import random 

13import os 

14 

15import numpy as np 

16 

17import gamdpy as gp 

18 

19# Some system parameters 

20nx, ny, nz = 6, 6, 10 

21rhoWall = 1.0 

22rhoFluid = 0.7 

23 

24# Setup a default fcc configuration 

25configuration = gp.Configuration(D=3) 

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

27configuration['m'] = 1.0 

28 

29 

30# Fluid particles have type '0', wall particles '1', dummy particles '2' 

31nwall, npart = 0, configuration.N 

32hlz = 0.5 * configuration.simbox.get_lengths()[2] 

33for n in range(npart): 

34 if configuration['r'][n][2] + hlz < 3.0: 

35 configuration.ptype[n] = 1 

36 nwall = nwall + 1 

37 

38nfluid = np.sum(configuration.ptype == 0) 

39nfluidWanted = int(nfluid * rhoFluid / rhoWall) 

40 

41while nfluid > nfluidWanted: 

42 idx = random.randint(0, npart - 1) 

43 if configuration.ptype[idx] == 0: 

44 configuration.ptype[idx] = 2 

45 nfluid = nfluid - 1 

46 

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

48 

49# Tether specifications.  

50tether = gp.Tether() 

51tether.set_anchor_points_from_types(particle_types=[1], spring_constants=[300.0], configuration=configuration) 

52 

53# Add gravity force  

54grav = gp.Gravity() 

55grav.set_gravity_from_types(particle_types=[0], forces=[0.01], configuration=configuration) 

56 

57# Temp relaxation for wall particles 

58relax = gp.Relaxtemp() 

59relax.set_relaxation_from_types(particle_types=[1], temperature=[2.], 

60 relax_times=[0.01],configuration=configuration); 

61 

62# Set the pair interactions 

63pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

64sig = [[1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, 0.0]] 

65eps = [[1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [0.0, 0.0, 0.0]] 

66cut = np.array(sig) * 2.5 

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

68 

69# Temperature 

70configuration.randomize_velocities(temperature=2.0) 

71 

72# Setup integrator: NVT 

73integrator = gp.integrators.NVE(dt=0.005) 

74 

75# Compute plan 

76compute_plan = gp.get_default_compute_plan(configuration) 

77 

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

79runtime_actions = [gp.TrajectorySaver(), 

80 gp.ScalarSaver()] 

81 

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

83sim = gp.Simulation(configuration, [pair_pot, tether, grav, relax], integrator, runtime_actions, 

84 num_timeblocks=100, steps_per_timeblock=64, 

85 storage='memory', compute_plan=compute_plan) 

86 

87prof = gp.CalculatorHydrodynamicProfile(configuration, 0) 

88 

89# Run simulation one block at a time 

90for block in sim.run_timeblocks(): 

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

92 prof.update() 

93 

94print(sim.summary()) 

95 

96prof.read() 

97 

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

99