From d0ec9520f238ef728ec7cfc53901f490cee08a1b Mon Sep 17 00:00:00 2001 From: Rafael Ravedutti Date: Thu, 24 Feb 2022 02:36:17 +0100 Subject: [PATCH] Write function to read PDB files and include data for Argon simulation Signed-off-by: Rafael Ravedutti --- data/anneal.mdp | 233 ++++++++++++++++++++++++++++++++++++++++ data/argon.top | 22 ++++ data/argon_start.pdb | 106 ++++++++++++++++++ data/heatup.mdp | 233 ++++++++++++++++++++++++++++++++++++++++ data/md.mdp | 233 ++++++++++++++++++++++++++++++++++++++++ gromacs/atom.c | 106 +++++++++++++++++- gromacs/includes/atom.h | 2 + 7 files changed, 933 insertions(+), 2 deletions(-) create mode 100644 data/anneal.mdp create mode 100644 data/argon.top create mode 100644 data/argon_start.pdb create mode 100644 data/heatup.mdp create mode 100644 data/md.mdp diff --git a/data/anneal.mdp b/data/anneal.mdp new file mode 100644 index 0000000..782f9d0 --- /dev/null +++ b/data/anneal.mdp @@ -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 diff --git a/data/argon.top b/data/argon.top new file mode 100644 index 0000000..2977bb1 --- /dev/null +++ b/data/argon.top @@ -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 diff --git a/data/argon_start.pdb b/data/argon_start.pdb new file mode 100644 index 0000000..67c4263 --- /dev/null +++ b/data/argon_start.pdb @@ -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 diff --git a/data/heatup.mdp b/data/heatup.mdp new file mode 100644 index 0000000..65c8e23 --- /dev/null +++ b/data/heatup.mdp @@ -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 diff --git a/data/md.mdp b/data/md.mdp new file mode 100644 index 0000000..e0dddb1 --- /dev/null +++ b/data/md.mdp @@ -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 diff --git a/gromacs/atom.c b/gromacs/atom.c index 74afc7a..907297a 100644 --- a/gromacs/atom.c +++ b/gromacs/atom.c @@ -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(¶m->input_file[len - 4], ".pdb", 4) == 0) { return readAtom_pdb(atom, param); } + if(strncmp(¶m->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; } diff --git a/gromacs/includes/atom.h b/gromacs/includes/atom.h index 22c5c5f..70dc7af 100644 --- a/gromacs/includes/atom.h +++ b/gromacs/includes/atom.h @@ -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*);