Hello,<br>
<br>
Thanks for your reply. As for the simulation box, I did not change anything and I am still saying that I am having a cubic box of size 30nm. When I view the molecule in ngmx I see the cubic box.. no problem with that. What I was referring to was the option on menu bar where one can select type of box to view the trajecotry and using say triclinic box I see chain fits better (almost completely)in the box while in the case of cubic (which is my real box) I see molecules is not fitted completely inside the box.. I was just asking why is this happening. as you said "none of that probably means anything" . I think also it is unnecessary for now so please dont get confused. <br>
<br>As you requested I am including ALL the commands I used from the very beginning starting from replicating a 4C unit. <br>
<pre style="text-align: justify;">editconf -f Eth4C.pdb -o Eth4C-d.gro -princ -d 0.058 -bt cubic
editconf -f Eth4C-d.gro -o Eth4C-d-index.gro -n index.ndx -princ<br>genconf -f Eth4C-d-index.gro -o Eth4C-d-index-rep30.gro -nbox 30 1 1 now I have PE chain of 60 units.<br><br>using proper rtp and hdb files:
pdb2gmx -f
Eth4C-d-index-rep30-res.gro -o Eth4C-d-index-rep30-res-pdb2.gro -p
Eth4C-d-index-rep30-res-pdb2.top -ff oplsaa >& output.pdb2gmx<p style="margin-bottom: 0cm;" align="JUSTIFY">Then I EM this last gro file but since dimensions are 15 * 0.5 0.5 nm I am getting the error message : </p>
The cut-off length is
longer than half the shortest box vector or longer than the smallest
box diagonal element. Increase the box size or decrease rlist.
<br>then increased the size >
<pre style="text-align: justify;">editconf -f Eth4C-d-index-rep30-res-pdb2.gro -o Eth4C-d-index-rep30-res-pdb2-boxsize30.gro -box 30 30 30 -angles 90 90 90<br><br>Then I EM the molecule with new box size 30nm and call this energy minimized structure PE60.gro. Now I want to replicate this chain. density demands 1 molecule in 5.76 nm3. As you said I forgot the density for a moment and took a large box.<br>
PE60 is already in a large enough box (30nm). this is all I did before doing any NPT.<br><br>As for compressing the system to get more realistic structure and approaching the density I am after, I did several NPT runs. All show system is not equilbrated and they all crash and box does not become smaller than 2.2 nm.<br>
<br>As you proposed I tried doing runs in stepwise manner. Below you see my trails:<br><br>First trial :P 120 bar for 100 ps then P80 bar for 500 ps then P50 for 500 ps. this last run crashes at 2.2 nm size<br><br>As Mark suggested I tried started with lower pressure:<br>
<br>Second trail: P 50 for 500ps > P40 bar for 500 ps > P35 for 1000ps > P30 for 1000 ps . Below is the results of this last run. before I had negative average Pressure but below the pressure is 30bar and real values fluctuate between -500 to 500 bar. <br>
T thermostat is working fine but LJ Columb terms not reaching equilibrium. Also box size is 2.7 nm which is larger than former trials.Molecule in not inside the box for almost the entire simulation.Only a small portion in inside and since it is long can not fit in.<br>
<br>Do i need to start with lower pressures and go to higher values. Because the size is not changing that much with 30 bar?<br><br>Energy Average RMSD Fluct. Drift Tot-Drift<br>-------------------------------------------------------------------------------<br>
Angle 736.582 30.0543 30.0543 -9.65991e-05 -0.0965993<br>Ryckaert-Bell. 147.603 13.5226 13.5185 -0.00115216 -1.15216<br>LJ-14 153.926 6.02125 6.02079 0.000258266 0.258266<br>
Coulomb-14 -106.141 1.90986 1.90613 0.00041316 0.413161<br>LJ (SR) -354.892 15.5047 15.4964 -0.00176674 -1.76675<br>Coulomb (SR) 217.529 2.63621 2.6302 -0.000615852 -0.615853<br>
Potential 794.607 37.4862 37.4764 -0.00295993 -2.95994<br>Kinetic En. 899.88 30.6454 30.6399 0.00201524 2.01525<br>Total Energy 1694.49 18.7071 18.7051 -0.000944686 -0.944688<br>
Temperature 299.806 10.2099 10.2081 0.000671404 0.671405<br>Pressure (bar) 30.0416 269.98 269.948 0.0145652 14.5652<br>Cons. rmsd () 5.11735e-06 3.49122e-07 3.49102e-07 1.29073e-11 1.29073e-08<br>
Box-X 2.79128 0.00907586 0.00906645 1.43079e-06 0.00143079<br>Box-Y 2.79128 0.00907586 0.00906645 1.43079e-06 0.00143079<br>Box-Z 2.79128 0.00907586 0.00906645 1.43079e-06 0.00143079<br>
Density (SI) 128.686 1.2601 1.25879 -0.000199256 -0.199256<br>pV 39.3202 353.469 353.428 0.0188014 18.8014<br>Vir-XX 327.606 298.203 298.201 -0.00356238 -3.56239<br>
Mu-X -0.028926 0.304356 0.30434 -1.08836e-05 -0.0108837<br>Lamb-System 1.00002 0.000368363 0.000368248 -3.18082e-08 -3.18083e-05<br>Heat Capacity Cv: 12.4935 J/mol K (factor = 0.00115974)<br>
<br>***********************************************************************************<br><br>itle = PE-Hexane <br>;define = -DPOSRES tells gromacs to perform position restrained dynamics/include posre.itp into topology used for position restraint<br>
pbc = xyz use priodic BCs in all directions<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 = 500000 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 = 100<br>nstfout = 0<br>nstlog = 100 frequency to write energies to log file<br>
nstxtcout = 10 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 ; tells gromacs how to model electrostatics. 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 where to start switching the Coulomb potential <br>rvdw-switch = 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.0 distance for coulomb cut-off <br>rvdw = 1.0 distance for coulomb cut-off <br>
<br>; Temperature coupling <br>Tcoupl = 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 = berendsen Pressure coupling is not on<br>
Pcoupltype = isotropic means the box expands and contracts in all directions (x,y,z) in order to maintain the proper pressure<br>tau_p = 0.5 time constant for coupling in ps <br>compressibility = 4.5e-5 compressibility of solvent used in simulation in 1/bar<br>
ref_p = 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 = all-bonds sets the LINCS constraint for all bonds<br>
constraint-algorithm = lincs<br><br><br><br><br><br><br>
<br><br></pre><p></p></pre><br>