Coverage for examples/generic_molecules.py: 98%

1 

2 

3import numpy as np 

4import gamdpy as gp 

5 

6# Sim. params  

7rho, temperature = 1.0, 1.5 

8NVE = False # If True -> k small 

9angle0, k = 2.0, 500.0 

10#rbcoef=[15.5000, 20.3050, -21.9170, -5.1150, 43.8340, -52.6070] 

11rbcoef=[.0, 50.0, .0, .0, .0, .0] 

12 

13# Generate configuration with a FCC lattice 

14configuration = gp.Configuration(D=3) 

15configuration.make_lattice(gp.unit_cells.FCC, cells=(8, 8, 8), rho=rho) 

16configuration['m'] = 1.0 

17configuration.randomize_velocities(temperature=temperature) 

18 

19# Make bonds 

20bond_potential = gp.harmonic_bond_function 

21bond_params = [[1.0, 1000.], ] 

22bond_indices = [] 

23for n in range(0, configuration.N, 4): 

24 bond_indices.append([n, n+1, 0]) 

25 bond_indices.append([n+1, n+2, 0]) 

26 bond_indices.append([n+2, n+3, 0]) 

27 

28bonds = gp.Bonds(bond_potential, bond_params, bond_indices) 

29 

30# Angles 

31angle_potential = gp.cos_angle_function 

32angle_params = [[k, angle0],] 

33angle_indices = [] 

34for n in range(0, configuration.N, 4): 

35 angle_indices.append([n, n+1, n+2, 0]) 

36 angle_indices.append([n+1, n+2, n+3, 0]) 

37 

38angles = gp.Angles(angle_potential, angle_indices, angle_params) 

39 

40# Dihedrals 

41dihedral_potential = gp.ryckbell_dihedral 

42dihedral_params = [rbcoef, ] 

43dihedral_indices = [] 

44for n in range(0, configuration.N, 4): 

45 dihedral_indices.append([n, n+1, n+2, n+3, 0]) 

46 

47dihedrals = gp.Dihedrals(dihedral_potential, dihedral_indices, dihedral_params) 

48 

49# Exlusion list 

50#exclusions = angles.get_exclusions(configuration) 

51#exclusions = bonds.get_exclusions(configuration) 

52exclusions = dihedrals.get_exclusions(configuration) 

53 

54 

55# Make pair potential 

56pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

57sig, eps, cut = 1.0, 1.0, 2.5 

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

59 

60# Make integrator 

61if NVE: 

62 integrator = gp.integrators.NVE(dt=0.001) 

63else: 

64 integrator = gp.integrators.NVT(temperature=temperature, tau=0.1, dt=0.002) 

65 

66# Compute plan 

67compute_plan = gp.get_default_compute_plan(configuration) 

68 

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

70runtime_actions = [gp.TrajectorySaver(), 

71 gp.ScalarSaver(), 

72 gp.MomentumReset(100)] 

73 

74#  

75 

76# Setup simulation 

77sim = gp.Simulation(configuration, [pair_pot, bonds, angles, dihedrals], integrator, runtime_actions, 

78 num_timeblocks=10, steps_per_timeblock=256, 

79 compute_plan=compute_plan, storage='memory') 

80 

81angles_array, dihedrals_array = [], [] 

82for block in sim.run_timeblocks(): 

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

84 angles_array.append( angles.get_angle(10, configuration) ) 

85 dihedrals_array.append( dihedrals.get_dihedral(10, configuration) ) 

86 

87print(sim.summary()) 

88 

89columns = ['U', 'W', 'K',] 

90data = np.array(gp.extract_scalars(sim.output, columns, first_block=1)) 

91temp = 2.0/3.0*np.mean(data[2])/configuration.N 

92Etot = data[0] + data[2] 

93Etot_mean = np.mean(Etot)/configuration.N 

94Etot_std = np.std(Etot)/configuration.N 

95 

96print("Temp: %.2f Etot: %.2e (%.2e)" % (temp, Etot_mean, Etot_std)) 

97print("Angle mean: %.2f (standard deviation %.2f) " % (np.mean(angles_array), np.std(angles_array))) 

98print("Dihedral mean: %.2f (standard deviation %.2f) " % (np.mean(dihedrals_array), np.std(dihedrals_array))) 

99