Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
#################################
#### PARAMETERS ########
#################################
#temperature of the system
variable Temp equal 1.0
variable tdamp equal 1.0
#Equilibration
variable run1 equal 10000
#Restart every these timesteps
variable trestart equal 10000
#Dump data every these timesteps
variable dumpfreq equal 100
#Dump trajectory every these timesteps
variable dumptraj equal 100
variable seedthermo equal 125485
################################
#### DEFINTIONS ########
################################
#Where the restart files of the run (non the equilibration) will be stored
variable folder index DUMPS
shell mkdir ${folder}
#The starting configuration is always the same. Also the name of the restart files
variable rname index Eq.5_1.N1000.data
variable simname index Run.R
units lj
atom_style angle
boundary p p p
neighbor 1.2 bin
neigh_modify every 1 delay 1 check yes
restart ${trestart} ${folder}/Restart.${simname}.
read_data ${rname}
####################################
#### ANGLE #######
####################################
angle_style cosine
angle_coeff 1 5.0 #20 sigma for realistic DNA (1 bead = 2.5 nm)
#######################################################
#### PAIRS -- REPULSIVE + TOPO II #######
#######################################################
pair_style hybrid lj/cut 1.12246 soft 1.12246
pair_modify shift yes
pair_coeff * * lj/cut 1.0 1.0 1.12246
## NOW THERE IS TOPO2 -- TYPE 4 BEADS ARE SOFT AND CAN BE CROSSED
pair_coeff * 4 soft 2.0 1.12246
pair_coeff * 5 soft 20.0 1.12246
###################################
#### FENE ##################
##################################
bond_style fene
#It turns off the LJ part of the FENE-bond potential between 1-2 atoms. However, the pairwise LJ interaction is still on
special_bonds fene#<=== I M P O R T A N T (new command).
bond_coeff 1 30.0 1.6 1.0 1.0
bond_coeff 2 30.0 1.6 1.0 1.0
####################################
#### THERMO AND INTEGRATION
####################################
thermo 1000
thermo_style custom step temp epair vol cpu
timestep 0.01
####################################
#### DUMPS #######
####################################
dump 2 all custom ${dumpfreq} ${folder}/${simname}.* id type mol mass x y z ix iy iz vx vy vz fx fy fz
dump_modify 2 sort id
##XYZ
dump 4 all xyz ${dumptraj} ${folder}/traj.xyz
dump_modify 4 element O N Ca P Au
compute MSDs all msd
compute GYRs all gyration
fix imprimir all print ${dumpfreq} "$(step) $(c_GYRs) $(c_MSDs[4])" screen no append radius_of_gyration.txt
#Compute total energies
compute ekin all ke
compute epot all pe
variable ekin equal c_ekin
variable epot equal c_epot
variable etot equal c_ekin+c_epot
fix energia all print ${dumpfreq} "$(step) ${ekin} ${epot} ${etot}" screen no append energy_total.txt
####################################
#### FIXES #######
####################################
fix 1 all nve
fix 2 all langevin ${Temp} ${Temp} ${tdamp} ${seedthermo}
run ${run1}