Regarding NPT ensemble

Dear All

I am trying to calculate 555 nm box of electrolyte solution. The solutes are mixed with solven by considering 1M concentration. However, to gain bulk properties I need to get equilibrium system at 1 atm pressure. Therefore, First I run the NVE with langevine fix. Then I use press/berendsen and temp/berendsen to get the 1 atm pressure. Then I run the npt. But I really got confuse with some issues

  1. When I run press/berendsen and temp/berendsen after NVE I got the my desired temp and pressure with a few steps. After that all values got constant ( e.g Pot_Eng, Tot_Eng everything). I am wondering is this OK for press/berendsen and temp/berendsen? ( Since I could not find any example for this fix I got confused)

  2. After press/berendsen and temp/berendsen calcualtion when I run npt with pressure 1 atm, my pressure changed with a large amount ( e.g from 1.02 to -895.23 ) and then pressure fluctuate ± 600 at every 100 steps.

*** Units is real** and command is like this

fix 1 nve all
fix 2 all langevin 303 303 100 67891

Then for berendsen

fix 1 all press/berendsen iso 1.0 1.0 1000 modulus 1000
fix 2 all temp/berendsen 303 303 100

Then for npt calculation

fix 1 all npt temp 303 303 100 iso 1.0 1.0 1000

I reach at this stage of calculation by using manual and mail lists. And I also tried to find a reasonable solution for present stage. However, I failed to make it, therefore, it will be very helpful if experts show me the way thus I could solve the problem or sort out my mistakes.

Thank in advance

Regard

Anderson

Houston University

I assume you unfixed all those commands between the successive runs.

Pressures fluctuate with npt or berendsen. If they are not fluctuating with

berendsen then something is wrong with your damping or dynamics.

What should occur is the time average should converge to your

dersired pressure.

Steve

Dear Steve

Yes, I unfixed those commands.

Pressure fluctuate with npt and for berendsen total constant. Confusions are began at this stage. I mentioned my input and output.

units real

boundary p p p

box tilt large

atom_style full

bond_style harmonic

angle_style harmonic

dihedral_style harmonic

pair_style lj/cut/coul/long 12.0 12.0

pair_modify mix geometric

special_bonds lj/coul 0.0 0.0 0.5

kspace_style pppm 0.0001

----------------- Atom Definition Section -----------------

read_data system.data

----------------- Settings Section -----------------

include “system.in.settings”

----------------- Run Section --------------------

group sodium type 1
group phosphorus type 7 8
group ringcarbon type 3
group ringoxygen type 4
group carbonyl_carbon type 5
group carbonyl_oxygen type 6
group hydrogen type 2

group EC type 2 3 4 5 6

write_data newdata.data

Dear Steve

Yes, I unfixed those commands.

Pressure fluctuate with npt and for berendsen total constant. Confusions are
began at this stage. I mentioned my input and output.

why are you setting the bulk modulus for fix press/berendsen to 1000?
this is equivalent to changing Pdamp to 10fs, which will significantly
suppress fluctuations.

axel.

Dear Axel

Thank you very much for your comments

If I use modulus less than 1000 then error comes like this

ERROR on proc 0: Non-numeric box dimensions - simulation unstable (…/pppm_cg.cpp:288)

Therefore, I used the bulk modulus higher value.

Further comments will be very helpful

Thank you very much

Regard

Jeams

This sounds like your simulation is unstable, unless you use an exceedingly high bulk modulus. A system that is close enough to equilibrium shouldn’t need that. Why don’t you try to relax the system first with a constant-volume ensemble, perhaps even without temperature coupling? If your initial configuration is unstable with NVE, there isn’t much point in trying other ensembles.

Giacomo

Hi Giacomo

Thank you very much for your comments

First I tried with NVE for 0.5 ns with timestep 1 and then go for */berendsend ensemble. However, phenomenon is similar like before. With modulus value below 1000 its shows same error massage. And same output for modulus value 1000.

Any further comments will be very helpful

Thanks in advance

Regard
Jeams