Dear Dr.Chaban,<br><br>I tried to get the same results as you got for density and total energy starting from the 8 chains is the box. I used the revised mdp file and since PME option takes much time I used shift instead for a moment to verify if I can get the density you obtained. Also I used constraints for bonds. density I get after 1800ps is 850 kg/m3 while what you got is a slightly different (after 1000ps density varies around 880!). Also the total energy you have is bout 2300 whereas I get somethong around 12000!. I dont know why they are so different..<br>
<br>In the second trial I used "constraints = none" just the "shift" for electrostatics is not the same as the mdp file you suggested. again results are not the same for total energy (20000 Kj/mol). Below is the average values...also simulation crashes at 1550 ps and the last density I am getting is 834! still less than what you got. dont really know why results are not the same. The reason simulation crashes could be due to dt=0.002 while removing constraints? <br>
<br>I have two more inquiries:<br><br>1- What I see is that the system is getting compressed to a size of around 3 3 3 nm is the first 20ps abruptly and in the remainder of time system reaches equilibrium... This is very different approach than what I was trying doing before. I used berendsen scheme for pressure coupling and box reduced is size slowly but I had difficulty reaching the size I want and most of the simulations crashed. Does this mean berenden is not suitable for compressing system for all systems?<br>
<br>2- You said I dont need NVT at all. But if I want to study the system I get from NPT giving the density I want, can I not do NVT for a fixed system size and calcualte energies or other quantities I am after?<br><br>Thanks for your attention<br>
Moeed<br><br><br>Energy Average RMSD Fluct. Drift Tot-Drift<br>-------------------------------------------------------------------------------<br>Bond 3885.23 141.365 139.725 0.0478427 74.3668<br>
Angle 5864.23 149.886 149.885 -0.00114897 -1.78597<br>Ryckaert-Bell. 2243.62 144.068 134.737 -0.113664 -176.679<br>LJ-14 1167.5 30.8742 30.6402 -0.00845667 -13.1451<br>
Coulomb-14 -553.954 31.3413 30.7865 -0.0130831 -20.3364<br>LJ (SR) -4898.73 405.434 386.584 -0.272304 -423.27<br>Coulomb (SR) 1250.22 56.2361 56.1252 0.0078672 12.2288<br>
Potential 8958.11 379.038 344.366 -0.352947 -548.621<br>Kinetic En. 10834.2 156.085 156.066 -0.00535835 -8.32903<br>Total Energy 19792.3 414.016 381.523 -0.358305 -556.95<br>
Temperature 300.069 4.32299 4.32248 -0.000148407 -0.230684<br>Pressure (bar) 29.5896 1263.17 1263.16 0.0112876 17.5454<br>Box-X 3.67408 1.72247 1.70558 -0.000536305 -0.833634<br>
Box-Y 3.44445 1.61482 1.59898 -0.000502786 -0.781532<br>Box-Z 2.45126 1.66735 1.64995 -0.000535423 -0.832263<br>Density (SI) 815.504 72.8491 70.7123 0.0390324 60.6721<br>
pV 54.9014 2111.1 2111.07 -0.0208174 -32.3586<br>Vir-XX 3606.98 1466.95 1466.92 0.0213665 33.2121<br>Box-Vel-XX -0.0182685 0.234428 0.232397 6.8633e-05 0.106683<br>
Mu-Y 0.102393 0.966426 0.966077 -5.79406e-05 -0.0900631<br>Heat Capacity Cv: 12.4757 J/mol K (factor = 0.000207551)<br><br><br>; Run control<br>integrator = md ; type of dynamics algorithm. Here md uses a leap-frog algorithm for integrating Newtons's eq of motion<br>
dt = 0.002 ; in ps !<br>nsteps = 900000 ; length of simulation= nsteps*dt <br>nstcomm = 1 ; frequency for center of mass motion removal <br>
<br>; Output control<br>nstenergy = 100 ; frequency to write energies to energy file. i.e., energies and other statistical data are stored every 10 steps<br>nstxout = 100 ; frequency to write coordinates/velocity/force to output trajectory file. how often snapshots are collected= nstxout*dt<br>
nstvout = 0<br>nstfout = 0<br>nstlog = 1000 ; frequency to write energies to log file<br>nstxtcout = 0 ; frequency to write coordinates to xtc trajectory<br>
<br>; Neighbor searching<br>nstlist = 10 ; frequency to update neighbor list. Neighborlist will be updated at least every 10 steps. Manual p80<br>ns_type = grid ; make a grid in the box and only check atoms in neighboring grid cells when constructing a new neighbor list every nstlist steps<br>
<br>; Electrostatics/VdW <br>coulombtype = Shift ;PME rev ; tells gromacs how to model electrostatics. Shift: Coulomb/LJ potential is decreased over the whole range and forces decay smoothly to zero between <br>
vdw-type = Shift ; rcoulomb-switch/rvw-switch & rcoulomb/rvdw<br>rcoulomb-switch = 0.9 ;0 ; where to start switching the Coulomb potential <br>rvdw-switch = 0.9 ;0 ; where to start switching the LJ potential <br>
<br>; Cut-offs<br>rlist = 1.1 ; in nm. Cut-off distance for short-range neighbor list<br>rcoulomb = 1.1 ;1.0 ; distance for coulomb cut-off <br>rvdw = 1.0 ; distance for coulomb cut-off <br>
<br>; Temperature coupling <br>Tcoupl = v-rescale ;berendsen<br>tc-grps = System ;HEX ; groups to couple to thermostat; Berendsen temperature coupling is on in these groups <br>
tau_t = 0.1 ;0.1 ; time constant for T coupling <br>ref_t = 300 ;300 ; reference T for coupling. When you alter the T, don't forget to change the gen_temp for velocity generation<br>
<br>; Pressure coupling<br>Pcoupl = Parrinello-Rahman ;berendsen<br>Pcoupltype = semiisotropic ;isotropic ; isotropic: means the box expands and contracts in all directions(x,y,z)in order to maintain the proper pressure; semiisotropic: isotropic in x & y directions<br>
tau_p = 1 1 ;0.5 ; time constant for coupling in ps <br>compressibility = 4.5e-5 4.5e-5 ; compressibility of solvent used in simulation in 1/bar<br>ref_p = 30.0 30.0 ; reference P for coupling in bar<br>
<br>; Velocity generation Generate velocites is on at 300 K. Manual p155<br>gen_vel = yes ; generate velocites according to Maxwell distribution at T: gen_temp with random gen seed gen_seed<br>
gen_temp = 300.0 ; T for Maxwell distribution <br>gen_seed = 173529 ; used to initialize random generator for random velocities<br><br>; Bonds<br>constraints = none ;all-bonds ; sets the LINCS constraint for all bonds<br>
constraint-algorithm = lincs<br><br>