initial velocities for MD after a normal mode calculation

[kai]

Dear Prof. Graeme and others,

i have a thought of using the results of wonderful normal mode calculation to run a molecular dynamics simulation. I plan to excite a specific mode of a molecule adsorbed on surface.

based on my understanding, it goes as follows:

1 getting an equilibrium adsorbed geometry
2 run dyn matrix to get the normal mode frequency
3 pick up the target normal mode with the vector for atom displacement
4 apply initial velocities to each atoms of the molecule in the equilibrium configuration
5 run MD

My question is,
how to find out the proper vector for the displacement of atoms from a certain mode of dyn matrix calculation. i thought the generated a movies could be something, but when i compared the two consecutive frames, i found the center of mass had moved....anyway of getting around?

[graeme]

The modes from the dynamical matrix code should be written to the modes.dat file. Displacing along the appropriate mode should be what you want. Alternatively, the dymmodes2xyz.pl could be used to generate a displacement along a particular mode. The lowest modes may correspond to displacement of the center of mass, if you have no frozen atoms, but presumably you want to displace along a higher mode. If exactly zeroing the CM is important, you can do this manually from the data in the modes.dat file.

[kai]

thanks, Graeme.

A question that sounds a bit off topic from this forum, but i am still assuming you are the right person to ask

In a MD calculation, would the overal momentum conserved during the course of the process? say, in the N, V, E conserved scheme (IBRION = 0, SMASS = -3)? I thought naturally this should come, but i tested a free CO molecule and to my surprise the overal momentum is not conserved ....how to understand this? do you have similar experience when desorbing CH4 ?

thx

[graeme]

In an NVE simulation, the momentum should be conserved. For NVT, it would depend upon the thermostat, but in general momentum would not be conserved.

With desorbing CH4, there are frozen atoms in the system and so there are the equivalent to external forces, and momentum is not conserved.

[kai]

thx, Graeme,

I thought about the fixed bulk atoms too in the NVE....but i still have the same problem even putting a single CO molecule....dont know where the problem is from...

a quick google guides me here, and you think that is relevant to vasp?

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2346486/

For a computer simulation the cost of evaluating forces can be very high. In particular, 2-body nonbonded interactions require O(N2) operations unless approximations are used. The most efficient approximations are obtained with methods such as particle–particle particle–mesh [1], particle–mesh Ewald (PME) [2], and multilevel summation [3–5], which split the potential into a short-range part calculated directly and a long-range part interpolated from a mesh. The use of a mesh yields an approximation Ua(R) to the long-range part of the potential energy which is not translation-invariant. This leads not only to a violation of Newton’s third law but causes particles to exert forces on themselves [6]!

[graeme]

Those references are not particularly relevant to DFT since they are based upon empirical potentials. There is a grid in vasp upon which the FFTs are done. A course grid will lead some some drift in the force, but this is generally small. You can check using LREAL=.FALSE. to reduce this drift, as well as increase the FFT grid (prec=accurate, or increase NG(X,Y,Z)X).

But this is all much less likely than more obvious problems. For example, is your time step small enough? Are you sure you are running NVE? Is your initial geometry reasonable, or are there initial high forces? Does the number of electronic iterations ever exceed your limit? ...