This cannot take place if you perform a constant-volume simulation, using a super cell of constant volume and imposing periodic boundary condition. When using a super cell, atoms will have the freedom to relax into their local equilibrium positions (the super cell volume and geometry can be held constant, while the coordinates of the atoms inside it are changed during MD steps). Also, you should check whether the Ewald sum (involved in the calculation of the ion-ion interaction energy) is calculated properly, since for long-range Coulomb interaction of the ionic point charges it is a conditionally-convergent series (i.e., inappropriately summed, it can diverge or oscillate between several limits on increasing the number of the summands). Other thing to consider, is the time-step used in the integration of the equation-of-motion; it should not be too large to get the system out of the Born-Oppenheimer energy surface.

Note added: Also, the K=0 Fourier component of potentials must be treated properly. By charge neutrality, these components corresponding to all potentials (except the exchange-correlation potential) can be set equal to zero (any other finite constant would do, since they cancel, reflecting charge neutrality).

Dear Behnam

I'm using the well-known pppm method for long range coulomb interactions. Also the time step is 0.1 fs. Surprisingly, the graphene oxide nanosheet stars to move and the volume increases constantly when the coulomb interactions included. additionally i use npt ensemble for structure relaxation and pbc was applied.

I just added a Note to my previous comment. Is charge neutrality properly accounted for in your simulations? Please see the added Note.

Also, please let know what basis set you are using. If not plane waves, there is always a problem arising from incompleteness of the basis set used, giving rise to the so-called Pulay forces, which are not physical.

How can i check if charge neutrality is properly accounted? sorry for stupid question i asked.

If you mean the total charge, it is -0.00001. There is too many atoms in my graphene oxide sheet so i could not use a high level quantum calculations to calculate the partial atomic charges so yes i used b3lyp/sto-3g and pop=esp of course gaussian 09.

As for charge neutrality, you have to check the Fourier components corresponding to K=0 (that is the unit-cell averages) of the external potential and the Hartree potential. The two must exactly cancel.

Note added: Also check the K=0 Fourier component of your number density. This (the average number of electrons per unit cell) must be exactly equal to the average number of positive charges in the unit cell.

All checked and is ok. I think, as you mentioned, Pulay force has a role. Which method you suggest for partial atomic charges calculations?

I would use plane waves, but that is a personal bias. There is also one other force that you may be neglecting. Unless the Kohn-Sham orbitals are self-consistent with respect to the underlying effective potential, the total force exerted on each atom is not equal to the Hellmann-Feynman force (the Hellmann-Feynman theorem applies to eignstates, and away from self-consistency the Kohn-Sham orbitals are not the eigenfunctions of the Kohn-Sham Hamiltonian). You can check for this by reducing the size of the MD steps (this is certainly relevant in such soft material as a graphene nanosheet which you are dealing with).

Thank you dear Behnam i'm trying to extract the atomic charges with different method hope to work.

You are welcome Sayyed Jalil. Did you find some problem with charge imbalance in your system?

No i didn't and i'm completely confused.

Dear Behnam is it technically right if i use nvt instead of npt? I wanna to extract the total vdw and coulomb interaction energy between 2, 3, ... layers of graphene oxide.

Regarding the first part of your question, see the following paper by Nosé (after whom the Nosé-Hoover thermostat has been named) and the references therein:

http://www.tandfonline.com/doi/abs/10.1080/00268978400101201

See also this:

http://williamhoover.info/Scans1960s/1967-1.pdf

Regarding the second part of your question, Van de Waals (vdW) forces cannot be accurately calculated within the framework of the local-density approximation, LDA; in this approximation, the long-distance behaviour of the exchange-correlation potential is described wrongly (for an isolated atom, it diminishes exponentially, and not correctly like a power law). There is an approximate formalism for compensating for this shortcoming, by Andersson, Langreth, and Lundqvist, which has been described here:

http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.76.102

See also:

http://physics.rutgers.edu/~langreth/reprints/ijqc05.pdf

For an application of this formalism, see:

http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.100.053002

Lastly, before simulating your intended system, I recommend you to resolve the problems of your computer program by simulating a simple system about which you know everything. This makes debugging of your program far more easier. I expect that the solution to your problem must be sought in the context of the details discussed by Nosé in the paper I have cited here above.

Dear Behnam, it seems some misunderstanding is here, i'm using LAMMPS code with classical MD simulation. I have just extracted the atomic partial charges using quantum mechanics calculations then classical MD simulation is done.

Dear Sayyed Jalil, I did not know the details of the code you are using (until now; now I know by briefly looking through the documentation of the code on the website of LAMMPS) and somehow your earlier comment indicating that your calculations are not quantum-mechanical had escaped my attention. I can therefore say nothing meaningful about your present problem. My earlier recommendation remains however, that to find out about the source of your problem it is advisable to consider a simple and transparent problem that you know everything about. For instance, take just a diatomic molecule and see how your simulation progresses, and if the atoms do not bind, why it is the case.

My problem unfortunately insists however thank you so much dear Behnam for your time and your attention.

With pleasure!

Hi Syyed, were you using NPT ensemble? Why not try NVE OR NVT to control the volume?

Hi dear Yuxiang, is it technically right if i use nvt instead of npt? I wanna to extract the total vdw and coulomb interaction energy between layers of graphene oxide.

Then the question becomes do you want to control the pressure in your system? If not, I believe NVT would be good. Since you were simulating two layers and with NPT, I assume that you don't want periodic boundary condition perpendicular to the plane, but you have to set the boundary as p p p in LAMMPS in order to use NPT. So the simulation box should be set much larger than the two layers.

If NPT drives your graphene sheets away, I guess there's something wrong with the coulomb parameter you were using, if so, there might be the same problem even with NVT.

Is it possible to do the simulation without periodic boundary? Try out and see what happens with and without coulombic interactions. The expansion may be due to coulombic repulsions.

Thank you dear Chin

I think you are right, i'll try that.