On Fri, Jul 11, 2014 at 10:27 PM, Justin Lemkul <jalem...@vt.edu> wrote:
> > > On 7/10/14, 9:30 PM, Rasoul Nasiri wrote: > >> On Fri, Jul 11, 2014 at 12:58 AM, Justin Lemkul <jalem...@vt.edu> wrote: >> >> >>> >>> On 7/10/14, 7:09 PM, Rasoul Nasiri wrote: >>> >>> Because >>>> >>>> 1- Estimation of T for different sub-systems say phases of gas, >>>> interface >>>> and liquid. >>>> Can be obtained with g_energy? >>>> >>>> >>>> No, not unless they are assigned as separate tc-grps in the .mdp file, >>> >> >> >> as sub-systems are changing during the simulation in terms of number of >> atoms in a non-equilibrium situation, I don't think g_energy can give >> relevant temperature values, even tc-grps are specified in mdp. Isn't? >> >> > No, you can't then. This also makes any sort of interpretation > non-trivial; if an atom dissociates from one group to which it was linked, > the interpretation of its new velocity when bonding to another group is now > somewhat discontinuous, is it not? > > > So it could make clrea why I did not select g_energy for determination of T and replied your question! > >> >> but then (1) you shouldn't be designing something affecting the physics >>> purely based on what you hope to analyze >>> >> >> >> >> >> >> and (2) you then already know the outcome. >>> >>> This is actually what g_traj can do; like I said, you just have to >>> account >>> for restraints. >>> >> >> >> can you clarify how it can be practically done? >> >> > You need to recalculate the number of degrees of freedom based on whatever > constraints there might be. You said there are none, so that may not be > the case. Please excuse the typo in the earlier message, "restraints" > should be "constraints." Which DOFs, do I need consider translational and rotational as well as vibrational ones? And how this number will affact crude values obtained by g_traj? Please consider constrained MD simulation done using OPLS. > > > >> >>> >>> 2- My trajectory is not constrain one as it has been recorded using >>> >>>> reactive FF. >>>> >>>> I need to estimate T with converted trajectory (reaxff---->gro). >>>> >>>> >>>> Converted trajectory? >>> >> >> >> yes, since non-reactive FF at high temperatures cannot be reliable due to >> vibrational effects, all my MD has been done by reactive one. >> >> >> Your first post suggested the run was done with Gromacs (being that you >>> had .tpr and .edr files). >>> >> >> >> I have produced .tpr file but not .edr. >> >> >> The outcomes you posted agree perfectly with constrained molecules, >>> >> >> >> The posted data refereed to un-reliability of g_traj for obtaining of T. I >> just wanted to know the reason via comparison with g_energy. Those >> obtained >> by non-reactive FF of OPLS. >> >> > This is where I was getting confused. Your very first message said you > had the "same trajectory and input files" but clearly this is not true. > You're comparing apples and oranges; an accurate statement of what's > really going on is very important and saves everyone time. > > That means all posted values in my first message obtained by ONLY OPLS FF just to understand how much g_traj and g_energy give us different T values on same trajectory and input file. > > >> but without more details on what you're actually doing, I'm not going to >>> hazard a guess. Recall that even if you set "constraints = none," SETTLE >>> is still applied to water molecules unless you use -DFLEXIBLE (which you >>> shouldn't for most normal purposes). >>> >>> >>> >>> My molecules are not water and g_traj gave following T values for same >> system using two different FFs >> >> 1- ReaxFF >> 0 399.321 >> 1 400.25 >> 2 401.08 >> 3 401.271 >> 4 400.67 >> 5 400.264 >> 6 400.752 >> 7 401.178 >> 8 401.391 >> 9 399.065 >> 10 401.78 >> >> As i said already I just converted reaxff trajectory to gro one and also >> used tpr file (constrain=none in mdp) for estimation of T with g_traj. >> >> 2- OPLS >> 0 269.581 >> 1 271.064 >> 2 271.309 >> 3 271.205 >> 4 271.499 >> 5 270.868 >> 6 272.247 >> 7 271.59 >> 8 271.614 >> 9 270.598 >> 10 271.662 >> >> >> How constrained MD results (OPLS) can be corrected? >> >> > Now you're losing me again. You're talking about converting some ReaxFF > trajectory without constraints, but were there constraints in the OPLS > trajectory? here there are two different values of T obtained by ReaxFF and OPLS. As you told there are no constrains in ReaxFF trajectory BUT there should be constrain in OPLS trajectory. > Can you give an actual listing of commands and .mdp files (if applicable) > for what you're doing? This is all very confusing. > > > Please just tell me how values of obtained T by g_traj can be corrected when I use trajectory and tpr files of constrain MD simulations (OPLS) as an input? command: g_traj_d -f md.trr -s topol -n n.ndx -ot T.xvg mdp file: ; RUN CONTROL PARAMETERS integrator = md ; Start time and timestep in ps tinit = 0 dt = 0.002 nsteps = 10000000 ;20ns ; For exact run continuation or redoing part of a run init_step = 0 ; Part index is updated automatically on checkpointing (keeps files separate) simulation_part = 1 ; mode for center of mass motion removal comm-mode = Linear ; number of steps for center of mass motion removal nstcomm = 10 ; group(s) for center of mass motion removal comm-grps = ; LANGEVIN DYNAMICS OPTIONS ; 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 (ps^2) for minimization of flexible constraints fcstep = 0 ; Frequency of steepest descents steps when doing CG nstcgsteep = 1000 nbfgscorr = 10 ; TEST PARTICLE INSERTION OPTIONS rtpi = 0.05 ; OUTPUT CONTROL OPTIONS ; Output frequency for coords (x), velocities (v) and forces (f) nstxout = 100000 ; 200ps nstvout = 100000 nstfout = 0 ; Output frequency for energies to log file and energy file nstlog = 5000 ;10ps nstcalcenergy = 10 nstenergy = 5000 ;10ps ; Output frequency and precision for xtc file nstxtcout = 5000 ;10ps xtc-precision = 5000 ; This selects the subset of atoms for the xtc file. You can ; select multiple groups. By default all atoms will be written. xtc_grps = System ; Selection of energy groups energygrps = ; NEIGHBORSEARCHING PARAMETERS ; nblist update frequency nstlist = 10 ; ns algorithm (simple or grid) ns_type = grid ; Periodic boundary conditions: xyz, no, xy pbc = xyz periodic_molecules = no ; nblist cut-off rlist = 1.5 ; long-range cut-off for switched potentials rlistlong = 1.5 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = PME rcoulomb-switch = 0 rcoulomb = 1.5 ; Relative dielectric constant for the medium and the reaction field epsilon_r = 1 epsilon_rf = 1 ; Method for doing Van der Waals vdw-type = switch ; cut-off lengths rvdw-switch = 1.1 rvdw = 1.3 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; Extension of the potential lookup tables beyond the cut-off table-extension = 1 ; Seperate tables between energy group pairs energygrp_table = ; 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 ; OPTIONS FOR WEAK COUPLING ALGORITHMS ; Temperature coupling tcoupl = v-rescale nsttcouple = -1 nh-chain-length = 10 ; Groups to couple separately tc-grps = System ; Time constant (ps) and reference temperature (K) tau_t = 0.1 ref_t = 400 ; Pressure coupling Pcoupl = parrinello-rahman Pcoupltype = isotropic nstpcouple = -1 ; Time constant (ps), compressibility (1/bar) and reference P (bar) tau_p = 4.0 compressibility = 13.31e-5 ; cyclopentane ref_p = 1.0 ; Scaling of reference coordinates, No, All or COM refcoord_scaling = No ; Random seed for Andersen thermostat andersen_seed = 815131 ; OPTIONS FOR BONDS constraints = all-bonds ; Type of constraint algorithm constraint_algorithm = Lincs ; Do not constrain the start configuration continuation = yes ; 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 = 6 ; 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 --------------------------------------------------------------------------- Rasoul > -Justin > > -- > ================================================== > > Justin A. Lemkul, Ph.D. > Ruth L. Kirschstein NRSA Postdoctoral Fellow > > Department of Pharmaceutical Sciences > School of Pharmacy > Health Sciences Facility II, Room 601 > University of Maryland, Baltimore > 20 Penn St. > Baltimore, MD 21201 > > jalem...@outerbanks.umaryland.edu | (410) 706-7441 > http://mackerell.umaryland.edu/~jalemkul > > ================================================== > -- > Gromacs Users mailing list > > * Please search the archive at http://www.gromacs.org/ > Support/Mailing_Lists/GMX-Users_List before posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.