Let assume that I would like to create amorphopus polymer sample composed of C32H66 (dotriacontane) chains. I am trying to do it in the following way:
1. Create a single molecule (e.g. in Avogadro), save the yxz coordinates.
2. Use e.g. fortran code and create from the xyz file the pdb file. Preview of some middle lines of the file is below (I see that formating is lost):
HETATM 31 C31 DTCN 1 31.632 16.485 -14.222 1.00 0.00 C
HETATM 32 C32 DTCN 1 33.154 16.573 -14.167 1.00 0.00 C
HETATM 33 H011 DTCN 1 0.328 -0.862 0.798 1.00 0.00 H
HETATM 34 H012 DTCN 1 0.344 -0.416 -0.951 1.00 0.00 H
3. Copy to my working catalog and modify the oplsaa.ff force field folder (as I have chosen this potential), especially aminoacids.rtp and aminoacids.hdb files (the names probably depends on the Gromacs distributions, I am using version 4.5.3):
aminoacids.rtp
; C32H66 Dotriacontane
[ DTCN ]
[ atoms ]
C01 opls_135 -0.180 1
C02 opls_136 -0.120 2
C03 opls_136 -0.120 3
C04 opls_136 -0.120 4
C05 opls_136 -0.120 5
C06 opls_136 -0.120 6
C07 opls_136 -0.120 7
C08 opls_136 -0.120 8
C09 opls_136 -0.120 9
C10 opls_136 -0.120 10
C11 opls_136 -0.120 11
C12 opls_136 -0.120 12
C13 opls_136 -0.120 13
C14 opls_136 -0.120 14
C15 opls_136 -0.120 15
C16 opls_136 -0.120 16
C17 opls_136 -0.120 17
C18 opls_136 -0.120 18
C19 opls_136 -0.120 19
C20 opls_136 -0.120 20
C21 opls_136 -0.120 21
C22 opls_136 -0.120 22
C23 opls_136 -0.120 23
C24 opls_136 -0.120 24
C25 opls_136 -0.120 25
C26 opls_136 -0.120 26
C27 opls_136 -0.120 27
C28 opls_136 -0.120 28
C29 opls_136 -0.120 29
C30 opls_136 -0.120 30
C31 opls_136 -0.120 31
C32 opls_135 -0.180 32
H011 opls_156 0.060 1
H012 opls_156 0.060 1
H013 opls_156 0.060 1
H021 opls_156 0.060 2
H022 opls_156 0.060 2
H031 opls_156 0.060 3
H032 opls_156 0.060 3
H041 opls_156 0.060 4
H042 opls_156 0.060 4
H051 opls_156 0.060 5
H052 opls_156 0.060 5
H061 opls_156 0.060 6
H062 opls_156 0.060 6
H071 opls_156 0.060 7
H072 opls_156 0.060 7
H081 opls_156 0.060 8
H082 opls_156 0.060 8
H091 opls_156 0.060 9
H092 opls_156 0.060 9
H101 opls_156 0.060 10
H102 opls_156 0.060 10
H111 opls_156 0.060 11
H112 opls_156 0.060 11
H121 opls_156 0.060 12
H122 opls_156 0.060 12
H131 opls_156 0.060 13
H132 opls_156 0.060 13
H141 opls_156 0.060 14
H142 opls_156 0.060 14
H151 opls_156 0.060 15
H152 opls_156 0.060 15
H161 opls_156 0.060 16
H162 opls_156 0.060 16
H171 opls_156 0.060 17
H172 opls_156 0.060 17
H181 opls_156 0.060 18
H182 opls_156 0.060 18
H191 opls_156 0.060 19
H192 opls_156 0.060 19
H201 opls_156 0.060 20
H202 opls_156 0.060 20
H211 opls_156 0.060 21
H212 opls_156 0.060 21
H221 opls_156 0.060 22
H222 opls_156 0.060 22
H231 opls_156 0.060 23
H232 opls_156 0.060 23
H241 opls_156 0.060 24
H242 opls_156 0.060 24
H251 opls_156 0.060 25
H252 opls_156 0.060 25
H261 opls_156 0.060 26
H262 opls_156 0.060 26
H271 opls_156 0.060 27
H272 opls_156 0.060 27
H281 opls_156 0.060 28
H282 opls_156 0.060 28
H291 opls_156 0.060 29
H292 opls_156 0.060 29
H301 opls_156 0.060 30
H302 opls_156 0.060 30
H311 opls_156 0.060 31
H312 opls_156 0.060 31
H321 opls_156 0.060 32
H322 opls_156 0.060 32
H323 opls_156 0.060 32
[ bonds ]
C01 C02
C02 C03
C03 C04
C04 C05
C05 C06
C06 C07
C07 C08
C08 C09
C09 C10
C10 C11
C11 C12
C12 C13
C13 C14
C14 C15
C15 C16
C16 C17
C17 C18
C18 C19
C19 C20
C20 C21
C21 C22
C22 C23
C23 C24
C24 C25
C25 C26
C26 C27
C27 C28
C28 C29
C29 C30
C30 C31
C31 C32
H011 C01
H012 C01
H013 C01
H021 C02
H022 C02
H031 C03
H032 C03
H041 C04
H042 C04
H051 C05
H052 C05
H061 C06
H062 C06
H071 C07
H072 C07
H081 C08
H082 C08
H091 C09
H092 C09
H101 C10
H102 C10
H111 C11
H112 C11
H121 C12
H122 C12
H131 C13
H132 C13
H141 C14
H142 C14
H151 C15
H152 C15
H161 C16
H162 C16
H171 C17
H172 C17
H181 C18
H182 C18
H191 C19
H192 C19
H201 C20
H202 C20
H211 C21
H212 C21
H221 C22
H222 C22
H231 C23
H232 C23
H241 C24
H242 C24
H251 C25
H252 C25
H261 C26
H262 C26
H271 C27
H272 C27
H281 C28
H282 C28
H291 C29
H292 C29
H301 C30
H302 C30
H311 C31
H312 C31
H321 C32
H322 C32
H323 C32
aminoacids.hdb
DTCN 32
3 4 H01 C01 C02 C03
2 4 H02 C02 C03 C04
2 4 H03 C03 C04 C05
2 4 H04 C04 C05 C06
2 4 H05 C05 C06 C07
2 4 H06 C06 C07 C08
2 4 H07 C07 C08 C09
2 4 H08 C08 C09 C10
2 4 H09 C09 C10 C11
2 4 H10 C10 C11 C12
2 4 H11 C11 C12 C13
2 4 H12 C12 C13 C14
2 4 H13 C13 C14 C15
2 4 H14 C14 C15 C16
2 4 H15 C15 C16 C17
2 4 H16 C16 C17 C18
2 4 H17 C17 C18 C19
2 4 H18 C18 C19 C20
2 4 H19 C19 C20 C21
2 4 H20 C20 C21 C22
2 4 H21 C21 C22 C23
2 4 H22 C22 C23 C24
2 4 H23 C23 C24 C25
2 4 H24 C24 C25 C26
2 4 H25 C25 C26 C27
2 4 H26 C26 C27 C28
2 4 H27 C27 C28 C29
2 4 H28 C28 C29 C30
2 4 H29 C29 C30 C31
2 4 H30 C30 C31 C32
2 4 H31 C31 C30 C29
3 4 H32 C32 C31 C30
4. Use pdb2gmx command to create .gro file:
pdb2gmx -f X.pdb
and select forcefield from the working catalog and no water
5. Use genbox command to create larger random sample:
genbox -ci conf.gro -nmol 30 -box 4.8 4.8 12 -o box.gro
6. Create (e.g. in Fortran) pdb file for this larger system:
HETATM 1 C01 DTCN 1 16.290 5.740 69.430 1.00 0.00 C
HETATM 2 C02 DTCN 1 17.270 6.150 68.320 1.00 0.00 C
HETATM 3 C03 DTCN 1 16.780 5.670 66.950 1.00 0.00 C
HETATM 4 C04 DTCN 1 17.770 6.080 65.850 1.00 0.00 C
HETATM 5 C05 DTCN 1 17.280 5.590 64.480 1.00 0.00 C
HETATM 6 C06 DTCN 1 18.250 6.000 63.370 1.00 0.00 C
HETATM 7 C07 DTCN 1 17.770 5.520 61.990 1.00 0.00 C
HETATM 8 C08 DTCN 1 18.760 5.930 60.890 1.00 0.00 C
HETATM 9 C09 DTCN 1 18.270 5.450 59.520 1.00 0.00 C
HETATM 10 C10 DTCN 1 19.250 5.860 58.420 1.00 0.00 C
HETATM 11 C11 DTCN 1 18.760 5.380 57.050 1.00 0.00 C
HETATM 12 C12 DTCN 1 19.740 5.790 55.940 1.00 0.00 C
HETATM 13 C13 DTCN 1 19.260 5.310 54.560 1.00 0.00 C
HETATM 14 C14 DTCN 1 20.230 5.720 53.460 1.00 0.00 C
HETATM 15 C15 DTCN 1 19.750 5.240 52.090 1.00 0.00 C
HETATM 16 C16 DTCN 1 20.720 5.650 50.980 1.00 0.00 C
HETATM 17 C17 DTCN 1 20.240 5.170 49.600 1.00 0.00 C
HETATM 18 C18 DTCN 1 21.220 5.590 48.500 1.00 0.00 C
HETATM 19 C19 DTCN 1 20.740 5.100 47.130 1.00 0.00 C
HETATM 20 C20 DTCN 1 21.710 5.520 46.020 1.00 0.00 C
HETATM 21 C21 DTCN 1 21.230 5.030 44.660 1.00 0.00 C
HETATM 22 C22 DTCN 1 22.200 5.450 43.550 1.00 0.00 C
HETATM 23 C23 DTCN 1 21.720 4.960 42.170 1.00 0.00 C
HETATM 24 C24 DTCN 1 22.700 5.380 41.070 1.00 0.00 C
HETATM 25 C25 DTCN 1 22.210 4.890 39.700 1.00 0.00 C
HETATM 26 C26 DTCN 1 23.190 5.310 38.590 1.00 0.00 C
HETATM 27 C27 DTCN 1 22.710 4.810 37.220 1.00 0.00 C
HETATM 28 C28 DTCN 1 23.690 5.230 36.110 1.00 0.00 C
HETATM 29 C29 DTCN 1 23.210 4.730 34.740 1.00 0.00 C
HETATM 30 C30 DTCN 1 24.180 5.160 33.630 1.00 0.00 C
HETATM 31 C31 DTCN 1 23.710 4.650 32.270 1.00 0.00 C
HETATM 32 C32 DTCN 1 24.670 5.070 31.170 1.00 0.00 C
HETATM 33 H011 DTCN 1 16.670 6.100 70.400 1.00 0.00 H
HETATM 34 H012 DTCN 1 16.200 4.640 69.480 1.00 0.00 H
HETATM 35 H013 DTCN 1 15.300 6.180 69.240 1.00 0.00 H
HETATM 36 H021 DTCN 1 18.270 5.720 68.550 1.00 0.00 H
HETATM 37 H022 DTCN 1 17.360 7.250 68.320 1.00 0.00 H
HETATM 38 H031 DTCN 1 15.780 6.110 66.740 1.00 0.00 H
HETATM 39 H032 DTCN 1 16.690 4.560 66.970 1.00 0.00 H
HETATM 40 H041 DTCN 1 18.770 5.640 66.060 1.00 0.00 H
HETATM 41 H042 DTCN 1 17.860 7.180 65.830 1.00 0.00 H
HETATM 42 H051 DTCN 1 16.280 6.030 64.260 1.00 0.00 H
HETATM 43 H052 DTCN 1 17.190 4.490 64.500 1.00 0.00 H
HETATM 44 H061 DTCN 1 19.260 5.570 63.580 1.00 0.00 H
HETATM 45 H062 DTCN 1 18.350 7.110 63.360 1.00 0.00 H
HETATM 46 H071 DTCN 1 16.770 5.960 61.790 1.00 0.00 H
HETATM 47 H072 DTCN 1 17.680 4.410 62.010 1.00 0.00 H
HETATM 48 H081 DTCN 1 19.750 5.490 61.110 1.00 0.00 H
HETATM 49 H082 DTCN 1 18.850 7.040 60.880 1.00 0.00 H
HETATM 50 H091 DTCN 1 17.270 5.890 59.310 1.00 0.00 H
HETATM 51 H092 DTCN 1 18.170 4.340 59.530 1.00 0.00 H
HETATM 52 H101 DTCN 1 20.250 5.420 58.620 1.00 0.00 H
HETATM 53 H102 DTCN 1 19.340 6.970 58.400 1.00 0.00 H
HETATM 54 H111 DTCN 1 17.760 5.820 56.830 1.00 0.00 H
HETATM 55 H112 DTCN 1 18.670 4.270 57.050 1.00 0.00 H
HETATM 56 H121 DTCN 1 20.740 5.350 56.150 1.00 0.00 H
HETATM 57 H122 DTCN 1 19.830 6.900 55.920 1.00 0.00 H
HETATM 58 H131 DTCN 1 18.250 5.750 54.350 1.00 0.00 H
HETATM 59 H132 DTCN 1 19.160 4.200 54.580 1.00 0.00 H
HETATM 60 H141 DTCN 1 21.240 5.280 53.660 1.00 0.00 H
HETATM 61 H142 DTCN 1 20.320 6.830 53.450 1.00 0.00 H
HETATM 62 H151 DTCN 1 18.740 5.680 51.880 1.00 0.00 H
HETATM 63 H152 DTCN 1 19.660 4.130 52.090 1.00 0.00 H
HETATM 64 H161 DTCN 1 21.730 5.220 51.190 1.00 0.00 H
HETATM 65 H162 DTCN 1 20.810 6.760 50.960 1.00 0.00 H
HETATM 66 H171 DTCN 1 19.240 5.610 49.390 1.00 0.00 H
HETATM 67 H172 DTCN 1 20.150 4.060 49.620 1.00 0.00 H
HETATM 68 H181 DTCN 1 22.210 5.150 48.710 1.00 0.00 H
HETATM 69 H182 DTCN 1 21.300 6.690 48.490 1.00 0.00 H
HETATM 70 H191 DTCN 1 19.740 5.540 46.920 1.00 0.00 H
HETATM 71 H192 DTCN 1 20.650 3.990 47.140 1.00 0.00 H
HETATM 72 H201 DTCN 1 22.710 5.090 46.230 1.00 0.00 H
HETATM 73 H202 DTCN 1 21.790 6.630 46.010 1.00 0.00 H
HETATM 74 H211 DTCN 1 20.230 5.470 44.450 1.00 0.00 H
HETATM 75 H212 DTCN 1 21.140 3.920 44.670 1.00 0.00 H
HETATM 76 H221 DTCN 1 23.210 5.020 43.750 1.00 0.00 H
HETATM 77 H222 DTCN 1 22.290 6.560 43.530 1.00 0.00 H
HETATM 78 H231 DTCN 1 20.720 5.390 41.960 1.00 0.00 H
HETATM 79 H232 DTCN 1 21.640 3.850 42.200 1.00 0.00 H
HETATM 80 H241 DTCN 1 23.700 4.950 41.280 1.00 0.00 H
HETATM 81 H242 DTCN 1 22.780 6.490 41.050 1.00 0.00 H
HETATM 82 H251 DTCN 1 21.210 5.320 39.480 1.00 0.00 H
HETATM 83 H252 DTCN 1 22.140 3.780 39.710 1.00 0.00 H
HETATM 84 H261 DTCN 1 24.190 4.880 38.800 1.00 0.00 H
HETATM 85 H262 DTCN 1 23.270 6.410 38.570 1.00 0.00 H
HETATM 86 H271 DTCN 1 21.710 5.230 37.010 1.00 0.00 H
HETATM 87 H272 DTCN 1 22.630 3.700 37.240 1.00 0.00 H
HETATM 88 H281 DTCN 1 24.700 4.810 36.330 1.00 0.00 H
HETATM 89 H282 DTCN 1 23.760 6.340 36.090 1.00 0.00 H
HETATM 90 H291 DTCN 1 22.210 5.150 34.530 1.00 0.00 H
HETATM 91 H292 DTCN 1 23.140 3.630 34.760 1.00 0.00 H
HETATM 92 H301 DTCN 1 25.190 4.740 33.850 1.00 0.00 H
HETATM 93 H302 DTCN 1 24.260 6.260 33.620 1.00 0.00 H
HETATM 94 H311 DTCN 1 22.700 5.080 32.040 1.00 0.00 H
HETATM 95 H312 DTCN 1 23.640 3.550 32.280 1.00 0.00 H
HETATM 96 H321 DTCN 1 24.750 6.180 31.110 1.00 0.00 H
HETATM 97 H322 DTCN 1 24.300 4.690 30.180 1.00 0.00 H
HETATM 98 H323 DTCN 1 25.680 4.640 31.340 1.00 0.00 H
HETATM 99 C01 DTCN 2 17.400 33.830 107.160 1.00 0.00 C
HETATM 100 C02 DTCN 2 18.030 32.720 106.310 1.00 0.00 C
...
7. Use pdb2gmx to obtain .gro and .topol files for this larger system
8. Prepare .mdp file for energy minimisation:
; LINES STARTING WITH ';' ARE COMMENTS
title = Minimization of Hen Egg White Lysozyme (1AKI.pdb) ADD YOUR NAME AND DATE HERE ; Title of run
; The following lines tell the program the standard locations where to find certain files
cpp = /lib/cpp ; Preprocessor
include = -I../top ; Directories to include in the topology format
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm options
; steep = steepest descent minimization
; MD = Leap Frog algorith for integrating Newtons equations of motion )
emtol = 10.0 ; Stop minimization when the energy changes by less than emtol kJ/mol.
nsteps = 1000 ; Maximum number of (minimization) steps to perform
nstenergy = 10 ; Write energies to disk every nstenergy steps
nstxtcout = 10 ; Write coordinates to disk every nstxtcout steps
xtc_grps = DTCN ; Which coordinate group(s) to write to disk
energygrps = DTCN ; Which energy group(s) to write to disk
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 5 ; Frequency to update the neighbor list and long range forces
ns_type = grid ; Method to determine neighbor list (simple, grid)
rlist = 1.4 ; Cut-off for making neighbor list (short range forces)
coulombtype = cut-off ; Treatment of long range electrostatic interactions
rcoulomb = 1.4 ; long range electrostatic cut-off
rvdw = 1.4 ; long range Van der Waals cut-off
constraints = none ; Bond types to replace by constraints
pbc = no ; Periodic Boundary Conditions (yes/no)
9. Modify .topol file:
1 opls_135 1 DTCN C01 1 -0.18 12.01
2 opls_136 1 DTCN C02 2 -0.12 12.01
3 opls_136 1 DTCN C03 3 -0.12 12.01
4 opls_136 1 DTCN C04 4 -0.12 12.01
5 opls_136 1 DTCN C05 5 -0.12 12.01
6 opls_136 1 DTCN C06 6 -0.12 12.01
7 opls_136 1 DTCN C07 7 -0.12 12.01
8 opls_136 1 DTCN C08 8 -0.12 12.01
9 opls_136 1 DTCN C09 9 -0.12 12.01
10 opls_136 1 DTCN C10 10 -0.12 12.01
11 opls_136 1 DTCN C11 11 -0.12 12.01
12 opls_136 1 DTCN C12 12 -0.12 12.01
13 opls_136 1 DTCN C13 13 -0.12 12.01
14 opls_136 1 DTCN C14 14 -0.12 12.01
15 opls_136 1 DTCN C15 15 -0.12 12.01
16 opls_136 1 DTCN C16 16 -0.12 12.01
17 opls_136 1 DTCN C17 17 -0.12 12.01
18 opls_136 1 DTCN C18 18 -0.12 12.01
19 opls_136 1 DTCN C19 19 -0.12 12.01
20 opls_136 1 DTCN C20 20 -0.12 12.01
21 opls_136 1 DTCN C21 21 -0.12 12.01
22 opls_136 1 DTCN C22 22 -0.12 12.01
23 opls_136 1 DTCN C23 23 -0.12 12.01
24 opls_136 1 DTCN C24 24 -0.12 12.01
25 opls_136 1 DTCN C25 25 -0.12 12.01
26 opls_136 1 DTCN C26 26 -0.12 12.01
27 opls_136 1 DTCN C27 27 -0.12 12.01
28 opls_136 1 DTCN C28 28 -0.12 12.01
29 opls_136 1 DTCN C29 29 -0.12 12.01
30 opls_136 1 DTCN C30 30 -0.12 12.01
31 opls_136 1 DTCN C31 31 -0.12 12.01
32 opls_135 1 DTCN C32 32 -0.18 12.01
33 opls_156 1 DTCN H011 1 0.06 1.01 (without modification I got here 33 instead of 1 in the 6th column)
34 opls_156 1 DTCN H012 1 0.06 1.01
35 opls_156 1 DTCN H013 1 0.06 1.01
36 opls_156 1 DTCN H021 2 0.06 1.01
37 opls_156 1 DTCN H022 2 0.06 1.01
10. Use grompp command to obtain tpr file:
grompp -f minim.mdp -c mol30.gro -p chg.top
11. Perform minimisation by:
mdrun
command
12. Copy confout.gro, chg.top and Jobscript file to another folder and prepare .mdp file for NVT relaxation:
; VARIOUS PREPROCESSING OPTIONS =
title = Yo
cpp = /lib/cpp
include = -I../top
define =
; RUN CONTROL PARAMETERS =
integrator = md
; start time and timestep in ps =
tinit = 0
dt = 0.002
nsteps = 25000
; number of steps for center of mass motion removal =
nstcomm = 1
comm-grps =
; LANGEVIN DYNAMICS OPTIONS =
; Temperature, friction coefficient (amu/ps) and random seed =
bd-fric = 0
ld-seed = 1993
; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout = 1000
nstvout = 0
nstfout = 0
; Output frequency for energies to log file and energy file =
nstlog = 1000
nstenergy = 1000
; Output frequency and precision for xtc file =
; This selects the subset of atoms for the xtc file. You can =
; select multiple groups. By default all atoms will be written. =
; Selection of energy groups =
energygrps = DTCN
; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist = 10
; ns algorithm (simple or grid) =
ns-type = grid
; Periodic boundary conditions: xyz or none =
pbc = xy
; nblist cut-off =
rlist = 1.4
; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype = Cut-off
rcoulomb = 1.4
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon-r = 1
; Method for doing Van der Waals =
vdw-type = Cut-off
; cut-off lengths =
rvdw-switch = 0
rvdw = 1.4
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = no
; Spacing for the PME/PPPM FFT grid =
fourierspacing = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used =
fourier-nx = 0
fourier-ny = 0
fourier-nz = 0
; EWALD/PME/PPPM parameters =
pme-order = 4
ewald-rtol = 1e-05
epsilon-surface = 0
; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling =
tcoupl = V-rescale
; Groups to couple separately =
tc-grps = System
; Time constant (ps) and reference temperature (K) =
tau-t = 0.1
ref-t = 300.0
; Pressure coupling =
Pcoupl = no
; SIMULATED ANNEALING CONTROL =
annealing = no
; Time at which temperature should be zero (ps) =
; GENERATE VELOCITIES FOR STARTUP RUN =
gen-vel = yes
gen-temp = 300
gen-seed = -1
; OPTIONS FOR BONDS =
constraints = all-bonds
; Type of constraint algorithm =
constraint-algorithm = Lincs
; Do not constrain the start configuration =
unconstrained-start = no
; Relative tolerance of shake =
shake-tol = 0.0001
; Highest order in the expansion of the constraint coupling matrix =
lincs-order = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs-warnangle = 30
; Convert harmonic bonds to morse potentials =
morse = no
; NMR refinement stuff =
; Distance restraints type: No, Simple or Ensemble =
disre = No
; Force weighting of pairs in one distance restraint: Equal or Conservative =
disre-weighting = Equal
; Use sqrt of the time averaged times the instantaneous violation =
disre-mixed = no
disre-fc = 1000
disre-tau = 0
; Output frequency for pair distances to energy file =
nstdisreout = 100
; Free energy control stuff =
free-energy = no
init-lambda = 0
delta-lambda = 0
sc-alpha = 0
sc-sigma = 0.3
; Non-equilibrium MD stuff =
acc-grps =
accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
energygrp-excl =
; Electric fields =
; Format is number of terms (int) and for all terms an amplitude (real) =
; and a phase angle (real) =
E-x =
E-xt =
E-y =
E-yt =
E-z =
E-zt =
; User defined thingies =
user1-grps =
user2-grps =
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0
13. Copy confout.gro, chg.top and Jobscript file to another folder and prepare .mdp file for NPT relaxation:
; VARIOUS PREPROCESSING OPTIONS =
title = Yo
cpp = /lib/cpp
include = -I../top
define =
; RUN CONTROL PARAMETERS =
integrator = md
; start time and timestep in ps =
tinit = 0
dt = 0.002
nsteps = 50000
; number of steps for center of mass motion removal =
nstcomm = 1
comm-grps =
; LANGEVIN DYNAMICS OPTIONS =
; Temperature, friction coefficient (amu/ps) and random seed =
bd-fric = 0
ld-seed = 1993
; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout = 1000
nstvout = 0
nstfout = 0
; Output frequency for energies to log file and energy file =
nstlog = 1000
nstenergy = 1000
; Output frequency and precision for xtc file =
; This selects the subset of atoms for the xtc file. You can =
; select multiple groups. By default all atoms will be written. =
; Selection of energy groups =
energygrps = DTCN
; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist = 10
; ns algorithm (simple or grid) =
ns-type = grid
; Periodic boundary conditions: xyz or none =
pbc = xy
; nblist cut-off =
rlist = 1.4
nwall = 2
wall-atomtype = C C
wall-type = 12-6
wall-r-linpot = -1.0
wall-ewald-zfac = 3
; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype = Cut-off
rcoulomb = 1.4
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon-r = 1
; Method for doing Van der Waals =
vdw-type = Cut-off
; cut-off lengths =
rvdw-switch = 0
rvdw = 1.4
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid =
fourierspacing = 0.16
; FFT grid size, when a value is 0 fourierspacing will be used =
fourier-nx = 0
fourier-ny = 0
fourier-nz = 0
; EWALD/PME/PPPM parameters =
pme-order = 4
ewald-rtol = 1e-05
epsilon-surface = 0
ewald-geometry = 3dc
; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling =
tcoupl = V-rescale
; Groups to couple separately =
tc-grps = System
; Time constant (ps) and reference temperature (K) =
tau-t = 0.1
ref-t = 300.0
; Pressure coupling =
Pcoupl = Parrinello-Rahman
Pcoupltype = Isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau-p = 2.0
compressibility = 4.5e-5
ref-p = 1.0
; SIMULATED ANNEALING CONTROL =
annealing = no
; Time at which temperature should be zero (ps) =
; GENERATE VELOCITIES FOR STARTUP RUN =
gen-vel = no
gen-temp = 300
gen-seed = -1
; OPTIONS FOR BONDS =
constraints = all-bonds
; Type of constraint algorithm =
constraint-algorithm = Lincs
; Do not constrain the start configuration =
unconstrained-start = yes
; Relative tolerance of shake =
shake-tol = 0.0001
; Highest order in the expansion of the constraint coupling matrix =
lincs-order = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs-warnangle = 30
; Convert harmonic bonds to morse potentials =
morse = no
; NMR refinement stuff =
; Distance restraints type: No, Simple or Ensemble =
disre = No
; Force weighting of pairs in one distance restraint: Equal or Conservative =
disre-weighting = Equal
; Use sqrt of the time averaged times the instantaneous violation =
disre-mixed = no
disre-fc = 1000
disre-tau = 0
; Output frequency for pair distances to energy file =
nstdisreout = 100
; Free energy control stuff =
free-energy = no
init-lambda = 0
delta-lambda = 0
sc-alpha = 0
sc-sigma = 0.3
; Non-equilibrium MD stuff =
acc-grps =
accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
energygrp-excl =
; Electric fields =
; Format is number of terms (int) and for all terms an amplitude (real) =
; and a phase angle (real) =
E-x =
E-xt =
E-y =
E-yt =
E-z =
E-zt =
; User defined thingies =
user1-grps =
user2-grps =
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0
Remarks:
A. I am using only pbc=xy periodic boundary conditions as I would like next make a bombardment of my sample, so I need surface and not broken (in the case of pbc) molecules in Z-direction.
B. Till the NPT relaxation this procedure works well. With this npt.mdp file the NPT relaxation also works, but I am not sure if such content of this file is the optimal one. I encountered a lot of error when trying some other options (e.g. PME algoritm): problem with spliting calclations to too many nodes (or threads), some segmentation fault errors, matrix inversion errors - all of them I found in the websites about Gromacs but I managed to run only with this configuration. The only error I have here is the "The domain decomposition grid has shifted too much in the Z-direction around cell 1 0 7" error (also reported somewhere in the web) but after restarting from checkpoint Gromacs managed to finish the calculations.
So, are my .mdp files (for minimisation, NVT and NPT relaxation) "ok" or something needs to be changed?
Thank you for reading all of this :)