H2 single molecule AIMD

[Lucas Lodeiro]

Hi all!

I am exploring single molecule AIMD (H2 in this case) for different purposes. The main problem that I face is temperature/rotational modes. I try to run a AIMD using NVE and NVT ensembles, but there are some problems.

For an NVE ensemble, the system starts rotating with a conserved angular moment, even if I set equal to zero the total angular moment... This disrupt the linear movement of H2 bond oscillation. The rotation couples with the bond and generates a fake two frequency oscillator. If, after the AIMD, I correct the trajectory to a zero total angular moment for all frames (simple for diatomic homonuclear molecule, but not so easy for polyatomic molecules), the oscillator recovers the single frequency motion.

For an NVT  ensemble, the thermalization is so hard for few atoms... and the system is very bad behaved.

Now I try an NVE ensemble with constrain MD, fixing the Y and Z position of both atoms, and put the bond along the X direction.

  &CONSTRAINT
    &FIXED_ATOMS
      COMPONENTS_TO_FIX      YZ
      LIST                   1 2
    &END FIXED_ATOMS
  &END CONSTRAINT

I think with it, I force the MD to evolve only over the X axis... but the initial step gives me a NaN temperature and kinetic energy, no matter what initial temperature I use. I attach both files.

The output file prints:

 Calculation of degrees of freedom
                                                      Number of atoms:         2
                                 Number of Intramolecular constraints:         4
                                 Number of Intermolecular constraints:         0
                                  Invariants(translation + rotations):         5
                                                   Degrees of freedom:        -3


 Restraints Information
                                  Number of Intramolecular restraints:         0
                                  Number of Intermolecular restraints:         0
 ************************** Velocities initialization **************************
 Initial Temperature                                                       NaN K
 COM velocity:                       NaN                 NaN                 NaN
 COM position:            5.669178398657      5.669178398657      5.669178398657
 Angular velocity:                   NaN                 NaN                 NaN
 *******************************************************************************

It seems, the degrees of freedom is: 3N - constrains - invariants...

Is there a way to force the MD only along one axis to avoid rotation?

Regards!

[Marcella Iannuzzi]

The translational and rotational degrees of freedom are subtracted, hence fixing 4 more degrees of freedom  by constraints you end up with a negative number of DoF, which is obviously not physical. 

Are you using PBC? The ANGVEL_ZERO keyword is ignored with PBC

Regards

Marcella

[Lucas Lodeiro]

Hi Marcella,

I try the MD with periodic and non periodic conditions. When ANGVEL_ZERO was applied for a non periodic cell, the simulation starts with zero angular velocity, but the numerical inaccuracies pumps some angular velocity at the first steps... and the molecule starts rotating.

As you mention, for a diatomic molecule it is not possible to apply any constraint... in all cases it will be not physical.

Is there a way to fix the rotational movement?

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/aed754ac-07de-40c2-9114-4370f7af7b16n%40googlegroups.com.

[Marcella Iannuzzi]

Hi Lucas, 

I would rather try to improve the accuracy of the integration of the equations of motion. 

The error in the forces is probably too large to guarantee the right behaviour of the conserved quantities. 

Regards

Marcella

[Lucas Lodeiro]

Hi Marcella an Travis,

I try both advice. To improve the forces, I increase the cutoff from 350 to 700 Ry and rel_cutoff from 40 to 60 Ry, and decrease the eps_scf (inner and outer scf) from 1.0E-6 to 1.0E-8. This does not solve the rotational problem, but the pump is lower than the no-improve input. The system shows equal the two fake peaks but the "distance" between them is smaller, and approaches the corrected trajectory with the true one peak.

For the second, I set angular velocity tolerance equal to 1.0E-11. In this case the rotation is off along all the trajectory, just a little temperature decreasing along the trajectory (17 K in 40 ps).

Firstly I want to thank you for your disposition and advice. Secondly, I have some questions.

In the case of improving the forces, as I am using a time step of 0.25 fs, I think the cutoff and eps values are the main parameters to improve the forces, isn't it? Do you expect a better improvement?

In the case of angular velocity tolerance, a value of 1.0E-11, is a too little one? There is a typical value to balance the angular movement and to subtract the lesser amount of kinetic energy?

Regards

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/df7eeb5f-ef2d-424c-a9c1-d8545a3012a7n%40googlegroups.com.