Write function to read PDB files and include data for Argon simulation

Signed-off-by: Rafael Ravedutti <rafaelravedutti@gmail.com>
This commit is contained in:
Rafael Ravedutti 2022-02-24 02:36:17 +01:00
parent ca7775a62a
commit d0ec9520f2
7 changed files with 933 additions and 2 deletions

233
data/anneal.mdp Normal file
View File

@ -0,0 +1,233 @@
;
; File 'mdout.mdp' was generated
; By user: bert (1001)
; On host: bertp3
; At date: Sat Dec 4 13:41:42 2004
;
; VARIOUS PREPROCESSING OPTIONS
include =
define =
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.002
nsteps = 500000
; For exact run continuation or redoing part of a run
init_step = 0
; mode for center of mass motion removal
comm-mode = Angular
; number of steps for center of mass motion removal
nstcomm = 1
; group(s) for center of mass motion removal
comm-grps =
; LANGEVIN DYNAMICS OPTIONS
; Temperature, friction coefficient (amu/ps) and random seed
bd-fric = 0
ld-seed = 1993
; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol = 10
emstep = 0.01
; Max number of iterations in relax_shells
niter = 20
; Step size (1/ps^2) for minimization of flexible constraints
fcstep = 0
; Frequency of steepest descents steps when doing CG
nstcgsteep = 1000
nbfgscorr = 10
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 5000
nstvout = 5000
nstfout = 0
; Checkpointing helps you continue after crashes
nstcheckpoint = 1000
; Output frequency for energies to log file and energy file
nstlog = 500
nstenergy = 500
; Output frequency and precision for xtc file
nstxtcout = 500
xtc-precision = 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =
; Selection of energy groups
energygrps = System
; long-range cut-off for switched potentials
rlistlong = -1
cutoff-scheme = Verlet
; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 50
; ns algorithm (simple or grid)
ns_type = grid
; Periodic boundary conditions: xyz (default), no (vacuum)
; or full (infinite systems only)
pbc = xyz
; nblist cut-off
rlist = 0.85
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = Cut-off
rcoulomb-switch = 0
rcoulomb = 0.85
; 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 = 0.85
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = Enerpres
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; 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
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no
; GENERALIZED BORN ELECTROSTATICS
; Algorithm for calculating Born radii
gb_algorithm = Still
; Frequency of calculating the Born radii inside rlist
nstgbradii = 1
; Cutoff for Born radii calculation; the contribution from atoms
; between rlist and rgbradii is updated every nstlist steps
rgbradii = 2
; Salt concentration in M for Generalized Born models
gb_saltconc = 0
; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
implicit_solvent = No
; 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 = 100
; Pressure coupling
Pcoupl = no
Pcoupltype = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; SIMULATED ANNEALING
; Type of annealing for each temperature group (no/single/periodic)
annealing = single
; Number of time points to use for specifying annealing in each group
annealing_npoints = 2
; List of times at the annealing points for each group
annealing_time = 0 1000
; Temp. at each annealing point, for each group.
annealing_temp = 100 25
; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel = yes
gen_temp = 100.0
gen_seed = 173529
; OPTIONS FOR BONDS
constraints = all-bonds
; Type of constraint algorithm
constraint-algorithm = Lincs
; Do not constrain the start configuration
unconstrained-start = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 1e-04
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter = 1
; 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
; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp_excl =
; NMR refinement stuff
; Distance restraints type: No, Simple or Ensemble
disre = No
; Force weighting of pairs in one distance restraint: Conservative or Equal
disre-weighting = Conservative
; 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
; Orientation restraints: No or Yes
orire = no
; Orientation restraints force constant and tau for time averaging
orire-fc = 0
orire-tau = 0
orire-fitgrp =
; Output frequency for trace(SD) to energy file
nstorireout = 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
; 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

22
data/argon.top Normal file
View File

@ -0,0 +1,22 @@
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 1 no 1.0 1.0
[ atomtypes ]
AR 39.948 0.0 A 0.00622127 9.69576e-06
[ moleculetype ]
; molname nrexcl
Ar 1
[ atoms ]
; id at type res nr residu name at name cg nr charge
1 AR 1 Ar Ar 1 0
[ system ]
; Name
Argon
[ molecules ]
; Compound #mols
AR 100

106
data/argon_start.pdb Normal file
View File

@ -0,0 +1,106 @@
HEADER Argon
REMARK THIS IS A SIMULATION BOX
CRYST1 77.395 77.395 77.395 90.00 90.00 90.00 P 1 1
MODEL 1
ATOM 1 Ar Ar 1 51.580 69.230 34.130 1.00 0.00
ATOM 2 Ar Ar 2 28.910 66.340 7.590 1.00 0.00
ATOM 3 Ar Ar 3 43.560 29.320 14.140 1.00 0.00
ATOM 4 Ar Ar 4 9.020 58.720 19.380 1.00 0.00
ATOM 5 Ar Ar 5 0.740 15.940 73.580 1.00 0.00
ATOM 6 Ar Ar 6 61.740 19.620 38.910 1.00 0.00
ATOM 7 Ar Ar 7 4.440 73.250 49.010 1.00 0.00
ATOM 8 Ar Ar 8 34.310 35.210 17.110 1.00 0.00
ATOM 9 Ar Ar 9 27.610 62.240 49.370 1.00 0.00
ATOM 10 Ar Ar 10 52.470 36.570 24.330 1.00 0.00
ATOM 11 Ar Ar 11 35.610 40.470 4.600 1.00 0.00
ATOM 12 Ar Ar 12 11.400 68.530 42.650 1.00 0.00
ATOM 13 Ar Ar 13 13.360 1.500 72.280 1.00 0.00
ATOM 14 Ar Ar 14 4.230 22.000 53.780 1.00 0.00
ATOM 15 Ar Ar 15 43.640 36.730 12.010 1.00 0.00
ATOM 16 Ar Ar 16 57.390 12.970 74.670 1.00 0.00
ATOM 17 Ar Ar 17 69.190 24.990 30.900 1.00 0.00
ATOM 18 Ar Ar 18 8.850 24.190 29.720 1.00 0.00
ATOM 19 Ar Ar 19 38.950 11.010 71.210 1.00 0.00
ATOM 20 Ar Ar 20 22.040 2.130 40.060 1.00 0.00
ATOM 21 Ar Ar 21 5.560 68.420 72.650 1.00 0.00
ATOM 22 Ar Ar 22 38.530 5.580 63.860 1.00 0.00
ATOM 23 Ar Ar 23 71.740 31.640 51.650 1.00 0.00
ATOM 24 Ar Ar 24 63.710 0.640 0.450 1.00 0.00
ATOM 25 Ar Ar 25 74.000 36.690 65.680 1.00 0.00
ATOM 26 Ar Ar 26 28.380 38.890 26.020 1.00 0.00
ATOM 27 Ar Ar 27 34.970 2.100 27.840 1.00 0.00
ATOM 28 Ar Ar 28 66.720 45.290 70.950 1.00 0.00
ATOM 29 Ar Ar 29 12.190 55.640 4.650 1.00 0.00
ATOM 30 Ar Ar 30 19.770 21.860 47.230 1.00 0.00
ATOM 31 Ar Ar 31 54.140 73.290 39.340 1.00 0.00
ATOM 32 Ar Ar 32 8.120 34.660 17.400 1.00 0.00
ATOM 33 Ar Ar 33 16.830 34.030 48.550 1.00 0.00
ATOM 34 Ar Ar 34 13.390 74.450 6.500 1.00 0.00
ATOM 35 Ar Ar 35 71.740 3.340 39.860 1.00 0.00
ATOM 36 Ar Ar 36 30.550 22.140 36.270 1.00 0.00
ATOM 37 Ar Ar 37 6.670 48.980 36.060 1.00 0.00
ATOM 38 Ar Ar 38 11.450 29.040 13.290 1.00 0.00
ATOM 39 Ar Ar 39 26.060 45.780 17.210 1.00 0.00
ATOM 40 Ar Ar 40 4.010 56.140 37.210 1.00 0.00
ATOM 41 Ar Ar 41 12.270 39.860 60.810 1.00 0.00
ATOM 42 Ar Ar 42 26.440 22.380 59.040 1.00 0.00
ATOM 43 Ar Ar 43 14.010 6.890 63.170 1.00 0.00
ATOM 44 Ar Ar 44 37.280 58.490 75.890 1.00 0.00
ATOM 45 Ar Ar 45 34.760 26.610 28.150 1.00 0.00
ATOM 46 Ar Ar 46 55.310 43.910 71.610 1.00 0.00
ATOM 47 Ar Ar 47 68.870 73.950 45.180 1.00 0.00
ATOM 48 Ar Ar 48 50.410 72.800 49.350 1.00 0.00
ATOM 49 Ar Ar 49 36.760 74.690 32.730 1.00 0.00
ATOM 50 Ar Ar 50 7.530 68.680 60.120 1.00 0.00
ATOM 51 Ar Ar 51 23.460 51.370 8.960 1.00 0.00
ATOM 52 Ar Ar 52 74.590 77.010 7.000 1.00 0.00
ATOM 53 Ar Ar 53 33.810 36.450 59.930 1.00 0.00
ATOM 54 Ar Ar 54 62.030 46.010 76.630 1.00 0.00
ATOM 55 Ar Ar 55 21.890 40.850 12.510 1.00 0.00
ATOM 56 Ar Ar 56 2.920 71.320 6.370 1.00 0.00
ATOM 57 Ar Ar 57 8.370 19.070 61.530 1.00 0.00
ATOM 58 Ar Ar 58 59.620 43.250 34.230 1.00 0.00
ATOM 59 Ar Ar 59 17.070 9.730 21.600 1.00 0.00
ATOM 60 Ar Ar 60 68.830 68.080 68.860 1.00 0.00
ATOM 61 Ar Ar 61 68.990 24.410 0.830 1.00 0.00
ATOM 62 Ar Ar 62 70.270 66.680 10.930 1.00 0.00
ATOM 63 Ar Ar 63 60.330 45.690 56.810 1.00 0.00
ATOM 64 Ar Ar 64 21.780 29.460 18.470 1.00 0.00
ATOM 65 Ar Ar 65 12.980 48.980 14.410 1.00 0.00
ATOM 66 Ar Ar 66 27.840 23.230 71.190 1.00 0.00
ATOM 67 Ar Ar 67 68.140 47.210 22.730 1.00 0.00
ATOM 68 Ar Ar 68 64.570 54.650 65.550 1.00 0.00
ATOM 69 Ar Ar 69 70.630 9.780 7.080 1.00 0.00
ATOM 70 Ar Ar 70 39.870 31.850 4.520 1.00 0.00
ATOM 71 Ar Ar 71 50.520 19.420 48.020 1.00 0.00
ATOM 72 Ar Ar 72 52.720 13.230 10.310 1.00 0.00
ATOM 73 Ar Ar 73 28.900 52.160 60.610 1.00 0.00
ATOM 74 Ar Ar 74 62.260 40.750 50.970 1.00 0.00
ATOM 75 Ar Ar 75 55.260 27.540 37.380 1.00 0.00
ATOM 76 Ar Ar 76 34.160 12.970 17.010 1.00 0.00
ATOM 77 Ar Ar 77 66.580 17.030 67.060 1.00 0.00
ATOM 78 Ar Ar 78 70.400 65.440 47.830 1.00 0.00
ATOM 79 Ar Ar 79 4.550 11.260 19.690 1.00 0.00
ATOM 80 Ar Ar 80 22.240 6.830 62.670 1.00 0.00
ATOM 81 Ar Ar 81 58.640 68.400 52.730 1.00 0.00
ATOM 82 Ar Ar 82 59.810 17.790 5.990 1.00 0.00
ATOM 83 Ar Ar 83 52.710 76.640 15.370 1.00 0.00
ATOM 84 Ar Ar 84 61.630 47.230 47.590 1.00 0.00
ATOM 85 Ar Ar 85 63.540 70.060 60.590 1.00 0.00
ATOM 86 Ar Ar 86 49.310 13.590 71.230 1.00 0.00
ATOM 87 Ar Ar 87 14.260 7.730 1.220 1.00 0.00
ATOM 88 Ar Ar 88 18.890 50.850 24.650 1.00 0.00
ATOM 89 Ar Ar 89 54.320 3.490 35.690 1.00 0.00
ATOM 90 Ar Ar 90 57.670 10.710 15.900 1.00 0.00
ATOM 91 Ar Ar 91 46.430 0.770 25.020 1.00 0.00
ATOM 92 Ar Ar 92 29.450 33.680 34.800 1.00 0.00
ATOM 93 Ar Ar 93 58.210 28.400 22.850 1.00 0.00
ATOM 94 Ar Ar 94 74.180 19.030 29.610 1.00 0.00
ATOM 95 Ar Ar 95 4.810 51.550 14.210 1.00 0.00
ATOM 96 Ar Ar 96 20.100 73.430 43.680 1.00 0.00
ATOM 97 Ar Ar 97 12.780 4.900 56.140 1.00 0.00
ATOM 98 Ar Ar 98 40.120 16.300 64.540 1.00 0.00
ATOM 99 Ar Ar 99 55.150 47.800 42.280 1.00 0.00
ATOM 100 Ar Ar 100 54.680 31.280 56.510 1.00 0.00
TER
ENDMDL

233
data/heatup.mdp Normal file
View File

@ -0,0 +1,233 @@
;
; File 'mdout.mdp' was generated
; By user: bert (1001)
; On host: bertp3
; At date: Sat Dec 4 13:41:42 2004
;
; VARIOUS PREPROCESSING OPTIONS
include =
define =
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.002
nsteps = 2500000
; For exact run continuation or redoing part of a run
init_step = 0
; mode for center of mass motion removal
comm-mode = Angular
; number of steps for center of mass motion removal
nstcomm = 1
; group(s) for center of mass motion removal
comm-grps =
; LANGEVIN DYNAMICS OPTIONS
; Temperature, friction coefficient (amu/ps) and random seed
bd-fric = 0
ld-seed = 1993
; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol = 10
emstep = 0.01
; Max number of iterations in relax_shells
niter = 20
; Step size (1/ps^2) for minimization of flexible constraints
fcstep = 0
; Frequency of steepest descents steps when doing CG
nstcgsteep = 1000
nbfgscorr = 10
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 5000
nstvout = 5000
nstfout = 0
; Checkpointing helps you continue after crashes
nstcheckpoint = 1000
; Output frequency for energies to log file and energy file
nstlog = 500
nstenergy = 500
; Output frequency and precision for xtc file
nstxtcout = 500
xtc-precision = 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =
; Selection of energy groups
energygrps = System
; long-range cut-off for switched potentials
rlistlong = -1
cutoff-scheme = Verlet
; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 50
; ns algorithm (simple or grid)
ns_type = grid
; Periodic boundary conditions: xyz (default), no (vacuum)
; or full (infinite systems only)
pbc = xyz
; nblist cut-off
rlist = 0.85
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = Cut-off
rcoulomb-switch = 0
rcoulomb = 0.85
; 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 = 0.85
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = Enerpres
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; 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
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no
; GENERALIZED BORN ELECTROSTATICS
; Algorithm for calculating Born radii
gb_algorithm = Still
; Frequency of calculating the Born radii inside rlist
nstgbradii = 1
; Cutoff for Born radii calculation; the contribution from atoms
; between rlist and rgbradii is updated every nstlist steps
rgbradii = 2
; Salt concentration in M for Generalized Born models
gb_saltconc = 0
; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
implicit_solvent = No
; 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 = 50
; Pressure coupling
Pcoupl = no
Pcoupltype = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; SIMULATED ANNEALING
; Type of annealing for each temperature group (no/single/periodic)
annealing = single
; Number of time points to use for specifying annealing in each group
annealing_npoints = 2
; List of times at the annealing points for each group
annealing_time = 0 5000
; Temp. at each annealing point, for each group.
annealing_temp = 50 125
; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel = no
gen_temp = 0.0
gen_seed = 173529
; OPTIONS FOR BONDS
constraints = all-bonds
; Type of constraint algorithm
constraint-algorithm = Lincs
; Do not constrain the start configuration
unconstrained-start = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 1e-04
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter = 1
; 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
; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp_excl =
; NMR refinement stuff
; Distance restraints type: No, Simple or Ensemble
disre = No
; Force weighting of pairs in one distance restraint: Conservative or Equal
disre-weighting = Conservative
; 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
; Orientation restraints: No or Yes
orire = no
; Orientation restraints force constant and tau for time averaging
orire-fc = 0
orire-tau = 0
orire-fitgrp =
; Output frequency for trace(SD) to energy file
nstorireout = 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
; 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

233
data/md.mdp Normal file
View File

@ -0,0 +1,233 @@
;
; File 'mdout.mdp' was generated
; By user: bert (1001)
; On host: bertp3
; At date: Sat Dec 4 13:41:42 2004
;
; VARIOUS PREPROCESSING OPTIONS
include =
define =
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.002
nsteps = 500000
; For exact run continuation or redoing part of a run
init_step = 0
; mode for center of mass motion removal
comm-mode = Linear
; number of steps for center of mass motion removal
nstcomm = 1
; group(s) for center of mass motion removal
comm-grps =
; LANGEVIN DYNAMICS OPTIONS
; Temperature, friction coefficient (amu/ps) and random seed
bd-fric = 0
ld-seed = 1993
; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol = 10
emstep = 0.01
; Max number of iterations in relax_shells
niter = 20
; Step size (1/ps^2) for minimization of flexible constraints
fcstep = 0
; Frequency of steepest descents steps when doing CG
nstcgsteep = 1000
nbfgscorr = 10
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 1000
nstvout = 1000
nstfout = 0
; Checkpointing helps you continue after crashes
nstcheckpoint = 1000
; Output frequency for energies to log file and energy file
nstlog = 100
nstenergy = 100
; Output frequency and precision for xtc file
nstxtcout = 500
xtc-precision = 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =
; Selection of energy groups
energygrps = System
; long-range cut-off for switched potentials
rlistlong = -1
cutoff-scheme = Verlet
; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 50
; ns algorithm (simple or grid)
ns_type = grid
; Periodic boundary conditions: xyz (default), no (vacuum)
; or full (infinite systems only)
pbc = xyz
; nblist cut-off
rlist = 0.85
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = Cut-off
rcoulomb-switch = 0
rcoulomb = 0.85
; 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 = 0.85
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = Enerpres
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; 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
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no
; GENERALIZED BORN ELECTROSTATICS
; Algorithm for calculating Born radii
gb_algorithm = Still
; Frequency of calculating the Born radii inside rlist
nstgbradii = 1
; Cutoff for Born radii calculation; the contribution from atoms
; between rlist and rgbradii is updated every nstlist steps
rgbradii = 2
; Salt concentration in M for Generalized Born models
gb_saltconc = 0
; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
implicit_solvent = No
; 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 = 100
; Pressure coupling
Pcoupl = no
Pcoupltype = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; SIMULATED ANNEALING
; Type of annealing for each temperature group (no/single/periodic)
annealing =
; Number of time points to use for specifying annealing in each group
annealing_npoints =
; List of times at the annealing points for each group
annealing_time =
; Temp. at each annealing point, for each group.
annealing_temp =
; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel = yes
gen_temp = 100.0
gen_seed = 173529
; OPTIONS FOR BONDS
constraints = all-bonds
; Type of constraint algorithm
constraint-algorithm = Lincs
; Do not constrain the start configuration
unconstrained-start = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 1e-04
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter = 1
; 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
; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
energygrp_excl =
; NMR refinement stuff
; Distance restraints type: No, Simple or Ensemble
disre = No
; Force weighting of pairs in one distance restraint: Conservative or Equal
disre-weighting = Conservative
; 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
; Orientation restraints: No or Yes
orire = no
; Orientation restraints force constant and tau for time averaging
orire-fc = 0
orire-tau = 0
orire-fitgrp =
; Output frequency for trace(SD) to energy file
nstorireout = 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
; 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

View File

@ -168,8 +168,109 @@ void createAtom(Atom *atom, Parameter *param)
}
}
int readAtom(Atom* atom, Parameter* param)
{
int type_str2int(const char *type) {
if(strncmp(type, "Ar", 2) == 0) { return 0; } // Argon
fprintf(stderr, "Invalid atom type: %s\n", type);
exit(-1);
return -1;
}
int readAtom(Atom* atom, Parameter* param) {
int len = strlen(param->input_file);
if(strncmp(&param->input_file[len - 4], ".pdb", 4) == 0) { return readAtom_pdb(atom, param); }
if(strncmp(&param->input_file[len - 4], ".dmp", 4) == 0) { return readAtom_dmp(atom, param); }
fprintf(stderr, "Invalid input file extension: %s\nValid choices are: pdb, dmp\n", param->input_file);
exit(-1);
return -1;
}
int readAtom_pdb(Atom* atom, Parameter* param) {
FILE *fp = fopen(param->input_file, "r");
char line[MAXLINE];
int read_atoms = 0;
if(!fp) {
fprintf(stderr, "Could not open input file: %s\n", param->input_file);
exit(-1);
return -1;
}
while(!feof(fp)) {
fgets(line, MAXLINE, fp);
char *item = strtok(line, " ");
if(strncmp(item, "CRYST1", 6) == 0) {
param->xlo = 0.0;
param->xhi = atof(strtok(NULL, " "));
param->ylo = 0.0;
param->yhi = atof(strtok(NULL, " "));
param->zlo = 0.0;
param->zhi = atof(strtok(NULL, " "));
param->xprd = param->xhi - param->xlo;
param->yprd = param->yhi - param->ylo;
param->zprd = param->zhi - param->zlo;
// alpha, beta, gamma, sGroup, z
} else if(strncmp(item, "ATOM", 4) == 0) {
char *label;
int atom_id, comp_id;
MD_FLOAT occupancy, charge;
atom_id = atoi(strtok(NULL, " ")) - 1;
while(atom_id + 1 >= atom->Nmax) {
growAtom(atom);
}
atom->type[atom_id] = type_str2int(strtok(NULL, " "));
label = strtok(NULL, " ");
comp_id = atoi(strtok(NULL, " "));
atom_x(atom_id) = atof(strtok(NULL, " "));
atom_y(atom_id) = atof(strtok(NULL, " "));
atom_z(atom_id) = atof(strtok(NULL, " "));
atom->vx[atom_id] = 0.0;
atom->vy[atom_id] = 0.0;
atom->vz[atom_id] = 0.0;
occupancy = atof(strtok(NULL, " "));
charge = atof(strtok(NULL, " "));
atom->ntypes = MAX(atom->type[atom_id] + 1, atom->ntypes);
atom->Natoms++;
atom->Nlocal++;
read_atoms++;
} else if(strncmp(item, "HEADER", 6) == 0 ||
strncmp(item, "REMARK", 6) == 0 ||
strncmp(item, "MODEL", 5) == 0 ||
strncmp(item, "TER", 3) == 0 ||
strncmp(item, "ENDMDL", 6) == 0) {
// Do nothing
} else {
fprintf(stderr, "Invalid item: %s\n", item);
exit(-1);
return -1;
}
}
if(!read_atoms) {
fprintf(stderr, "Input error: No atoms read!\n");
exit(-1);
return -1;
}
atom->epsilon = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
atom->sigma6 = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
atom->cutforcesq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
atom->cutneighsq = allocate(ALIGNMENT, atom->ntypes * atom->ntypes * sizeof(MD_FLOAT));
for(int i = 0; i < atom->ntypes * atom->ntypes; i++) {
atom->epsilon[i] = param->epsilon;
atom->sigma6[i] = param->sigma6;
atom->cutneighsq[i] = param->cutneigh * param->cutneigh;
atom->cutforcesq[i] = param->cutforce * param->cutforce;
}
fprintf(stdout, "Read %d atoms from %s\n", read_atoms, param->input_file);
fclose(fp);
return read_atoms;
}
int readAtom_dmp(Atom* atom, Parameter* param) {
FILE *fp = fopen(param->input_file, "r");
char line[MAXLINE];
int natoms = 0;
@ -258,6 +359,7 @@ int readAtom(Atom* atom, Parameter* param)
}
fprintf(stdout, "Read %d atoms from %s\n", natoms, param->input_file);
fclose(fp);
return natoms;
}

View File

@ -59,6 +59,8 @@ typedef struct {
extern void initAtom(Atom*);
extern void createAtom(Atom*, Parameter*);
extern int readAtom(Atom*, Parameter*);
extern int readAtom_pdb(Atom*, Parameter*);
extern int readAtom_dmp(Atom*, Parameter*);
extern void growAtom(Atom*);
extern void growClusters(Atom*);