A Rotational/Roto-translational Constraint Method for Condensed Matter
Jitai Yang Ke Li Jia Liu Jia Nie Hui Li

	
	
	
	
	
	
	
	
		
		
    < g r a p h i c s >

	
	
	
	
	
	
	





In condensed matter physics, particularly in perovskite materials, the rotational motion of molecules and ions  is associated with important issues such as ion conduction mechanism.
Constrained Molecular Dynamics (MD) simulations offer a means to separate translational, vibrational, and rotational motions, enabling the independent study of their effects. 
In this study, we introduce a rotational and roto-translational constraint algorithm based on  the Velocity Verlet integrator, which has been implemented into a homebrew version of the CP2K package. The MD results show that our program can selectively constrain the molecules and ions in the system and support long-time MD runs. The algorithm can help the future study of important rotation-related dynamic problems in condensed matter systems.
	
	
	
	
	


§ INTRODUCTION


Rotation of molecules/ions in condensed matter, especially in perovskite crystals, can have close relation with important problems: the molecular rotation can drive phase transition in amorphous ice<cit.>; the volumetric lattice strain and CH3NH3 molecular rotational degrees of freedom(DOF) can cooperate to create and stabilize polarons in perovskite solar cell materials<cit.>; in the study of solid electrolytes, anion rotation is associated with fast ion conduction, the coupling between anion rotation and cation diffusion is named as paddle-wheel mechanism which has drawn general interest<cit.>.

Experimental techniques can provide statistically averaged information on molecule/ion rotation. When combined with ab initio molecular dynamics(AIMD), further microscopic information can be obtained.

In addition, fixing the motion of molecules/ions in AIMD allows for an approximate calculation of the rotation effect on material properties<cit.>. However, in these studies, vibration and translation of the molecules/ions were also fixed. A more rational approach would be to fix only the rotation of the molecules/ions.

The Nola group derived a roto-translational constraint algorithm 

<cit.>. Their method has been implemented in GROMOS<cit.> and GROMACS<cit.>. 
They hope that the algorithm will improve the efficiency of molecular dynamics(MD) running and be used in molecular docking research. 

Their work is rigorous and original, however, the method is implemented from a biomolecule perspective and based on the leap frog integrator<cit.>. Most popular material MD software, including LAMMPS<cit.>, VASP<cit.>, and CP2K<cit.>, do not support the leap frog integrator, which hindered the application of the method to condensed matter. 

The application of rotational constraint in AIMD is also a problem. Most biomolecule simulation software is specialized in MD simulations using the empirical force field or simple QM/MM methods<cit.>, and do not support AIMD. The interfaces of biomolecule simulation software are improving, for instance, GROMACS can be combined with CP2K to use the quantum mechanical/molecular mechanical(QM/MM) method. However, it is still not as convenient and efficient as implementing the constraint algorithms directly in material simulation software. To apply the rotational constraint algorithm in AIMD, it is necessary to rely on software such as VASP, CP2K, etc.




The Krimm group derived a general constraint algorithm named WIGGLE<cit.>; an angular momentum constraint method was mentioned in their work.
The constrained accelerations are used in their derivation. They claimed that the WIGGLE algorithm is almost an order of magnitude more accurate than the RATTLE algorithm,  but the algorithm is rarely used, probably because it has not been implemented in popular dynamics software. 

To the best of our knowledge, there is no popular material simulation software that implements the rotational constraint algorithm. This work presents a rotational/roto-translational constraint algorithm based on the Velocity Verlet integrator<cit.>, which is commonly used in materials simulation software. The algorithm has been implemented into a homebrew version of the CP2K program, and has been tested in MD running. The results indicate that our program can selectively constrain the rotational/roto-translational motions of molecules in the system and support stable long-time runs.













§ THEORY




The theoretical derivation of the rotational constraints using the Velocity Verlet integrator is presented in this section. Our derivation follows the RATTLE algorithm<cit.>, which is widely applied in MD software. 
Constraints for linear molecules are special cases and are not discussed in this work.
Since the roto-translational constraints differ only in the initial conditions from the rotational constraints, their derivation is presented in the supporting information.



We describe the constrained atoms using the center-of-mass(COM) coordinates system for convenience:

    𝐫_a = 𝐑_a -∑_a=1^n_a m_a𝐑_a /∑_a=1^n_am_a,
       ∑_a=1^n_a m_a 𝐫_a=0.


n_a is the total number of constrained atoms, m_a is the mass of ath atom. 𝐑_a is the position vector of ath atom in Cartesian coordinates system, while 𝐫_a is the position vector of ath atom in the COM coordinates system.


We use the same rotational constraint conditions as Amadei  groups<cit.>.

    χ^(r) = ∑^n_a_a=1𝐫^0_a × m_a 𝐫_a(t) = 0.


The superscript r indicates the quantity corresponds to a positional constraint, and 𝐫_a^0 is the rotational reference position of ath atom.

Another option is to use condition equations without mass. 
We prefer the condition equations with mass because it is close to the rotational Eckart conditions which divide the nuclear motion into translations, rotations, and vibrations approximately<cit.>.


The constraints are applied to the whole molecule, some bond wiggling is allowed. In particular, if some atom's rotational reference position is at the COM(𝐫_a^0 = 0), the atom will be constrained in a manner that is related to the COM and has no direct contribution to the constraint conditions.

Differentiating the equations at time t yields the velocity constraints 

    χ^(v) = ∑^n_a_a=1𝐫^0_a × m_a 𝐯_a(t) = 0.

The superscript v indicates the quantity corresponds to a velocity constraint, 𝐯_a(t) is the velocity of ath atom at time t. The rotational reference positions are constants. If we replace 𝐫_a^0 with 𝐫_a(t), χ^(v) will be  the total angular moment of the constrained atoms.

Verlocity Verlet integrator has position part and velocity part

    𝐫_a(t+δ t)=𝐫_a(t)+δ t  𝐯_a(t)+1/2 m_aδ t^2 𝐟_a(t), 
       𝐯_a(t+δ t)=𝐯_a(t)+1/2 m_aδ t(𝐟_a(t)+𝐟_a(t+δ t)).

Consider the position part with the constraint force first

    𝐫_a(t+δ t)=𝐫_a^'(t+δ t)+1/2m_aδ t^2  𝐠_a^(r)(t).

𝐠_a^(r) is the positional constraint force on ath atom, the prime symbol in 𝐫_a^' indicates the quantities without constraints.

Multiply both sides of the equation by m_a, cross product both sides of the equation by 𝐫_a^0 from left, then sum over atoms, we get

    ∑^n_a_a=1𝐫_a^0 × m_a𝐫_a(t+δ t)=
    	∑^n_a_a=1𝐫_a^0 × 
    	[ m_a𝐫_a^'(t+δ t)+1/2δ t^2  𝐠_a^(r)(t)].

The left side of the equation is the same as our rotational  constraint conditions Eq. (<ref>), which should be zero at anytime. So

    ∑^n_a_a=1𝐫_a^0 × 
    	[ m_a𝐫_a^'(t+δ t)+1/2δ t^2  𝐠_a^(r)(t)] = 0.

The positional constraint force can be written as 

    𝐠_a^(r)=-∑_i=1^n_rλ_i^(r)∇_𝐫_a χ_i^(r).

n_r is the number of DOF are constrainted, n_r =3 in rotational constraint for nonlinear molecule.
λ_i^(r) is the ith Lagrange multipliers of position constraint. The minus symbol on the right side is just for definition because we can write it into the λ_i^(r). 

Substitute Eq. (<ref>) into Eq. (<ref>). Take (x_1, y_1, z_1, x_2, y_2, z_2, …, x_a,y_a,z_a) as atoms position in COM coordinate system, (X_1, Y_1, Z_1, X_2, Y_2, Z_2, …, X_a,Y_a,Z_a) as atoms position of rotational reference in COM coordinate system. We get

    g_a,x^(r) = -∑_i=1^n_rλ_i^(r)∂χ_i^(r)/∂ x_a
    		=m_a(- Z_aλ_2^(r)+ Y_aλ_3^(r)), 
       g_a,y^(r) = -∑_i=1^n_rλ_i^(r)∂χ_i^(r)/∂ y_a
    		=m_a(  Z_aλ_1^(r) -X_aλ_3^(r)), 
        g_a,z^(r) = -∑_i=1^n_rλ_i^(r)∂χ_i^(r)/∂ z_a
    		=m_a(-Y_aλ_1^(r) +X_aλ_2^(r)).

Define

    𝐜^(r)= ∑_a=1^n_a𝐫_a^0 × m_a𝐫_a^'(t+δ t).

Substitute Eq. (<ref>) and Eq. (<ref>) into Eq. (<ref>), we will get equations

    𝐜_x^(r) +∑_a=1^n_a1/2 m_aδ t^2
    		( -(Y_a^2 +Z_a^2) λ_1^(r)
    		+ X_aY_a λ_2^(r)
    		+ X_aZ_a λ_3^(r)) =0 
    𝐜_y^(r) +∑_a=1^n_a1/2 m_aδ t^2
    		(   X_aY_a λ_1^(r)
    		-(X_a^2 +Z_a^2) λ_2^(r)
    		+ Y_aZ_a λ_3^(r))  = 0 
    𝐜_z^(r) +∑_a=1^n_a1/2 m_aδ t^2
    		(    X_a Z_a λ_1^(r)
    		+Y_aZ_a λ_2^(r)
    		-(X_a^2 +Y_a^2)  λ_3^(r)) = 0.

In matrix form is

    𝐜^(r) + Θ^(r)λ^(r)= 0,
    
    		   λ^(r)= 
    		-(Θ^(r))^-1𝐜^(r).

Then we can solve all Lagrange multipliers λ^(r). From Eq. (<ref>) and Eq. (<ref>) we can get constraint force and constrained positions at the next time step 𝐫_a(t+δ t ).


The velocity part of Velocity Verlet integrator with constraint force is

    𝐯_a(t+δ t)=𝐯_a^'(t+δ t)+1/2m_aδ t  𝐠_a^(v)(t+δ t).

𝐠_a^(v) is the velocity constraint force on ath atom. As we have done in the position part, we multiply both sides of the equation by m_a, cross product both sides of the equation by 𝐫_a^0 from left, then sum over atoms

    ∑_a=1^n_a𝐫_a^0 × m_a 𝐯_a(t+δ t)
    	=∑_a=1^n_a𝐫_a^0 ×[m_a 𝐯_a^'(t+δ t)+1/2δ t 𝐠_a^(v)(t + δ t)].

From Eq. (<ref>)

    ∑_a=1^n_a𝐫_a^0 ×[ m_a
    	𝐯_a^'(t+δ t)+1/2δ t 𝐠_a^(v)(t+δ t)]
    	= 0.

Eq. (<ref>) is a semi-holonomic constraint condition.  The corresponding constraint force is<cit.>

    𝐠_a^(v) = - ∑ ^n_v _i=1λ _i^(v)∇_𝐯_a χ_i^(v)

Note that the gradient is with respect to the velocity 𝐯_a rather than positions. The rotational reference positions are constants, so the components of 𝐠_a^(v) has the same formalism as 𝐠_a^(r) in Eq. (<ref>)

    g_a,x^(v) = -∑_i=1^n_vλ_i^(v)∂χ_i^(v)/∂ v_a,x
    		=m_a(- Z_aλ_2^(v)+ Y_aλ_3^(v)), 
       g_a,y^(v) = -∑_i=1^n_vλ_i^(v)∂χ_i^(v)/∂ v_a,y
    		=m_a(  Z_aλ_1^(v) -X_aλ_3^(v)), 
        g_a,z^(v) = -∑_i=1^n_vλ_i^(v)∂χ_i^(v)/∂ v_a,z
    		=m_a(-Y_aλ_1^(v) +X_aλ_2^(v)).

Define

    𝐜^(v)= ∑_a=1^n_a𝐫_a^0 × m_a𝐯_a^'(t+δ t).

Substitute Eq. (<ref>) and Eq. (<ref>) into Eq. (<ref>) we will get 

    𝐜_x^(v) +∑_a=1^n_a1/2 m_aδ t
    		( -(Y_a^2 +Z_a^2) λ_1^(v)
    		+ X_aY_a λ_2^(v)
    		+ X_aZ_a λ_3^(v)) =0 
    𝐜_y^(v) +∑_a=1^n_a1/2 m_aδ t
    		(    X_aY_a λ_1^(v)
    		-(X_a^2 +Z_a^2) λ_2^(v)
    		+ Y_aZ_a λ_3^(v))  = 0 
    𝐜_z^(v) +∑_a=1^n_a1/2 m_aδ t
    		(    X_a Z_a λ_1^(v)
    		+Y_aZ_a λ_2^(v)
    		-(X_a^2 +Y_a^2)  λ_3^(v)) = 0.

In matrix form is

    𝐜^(v) + Θ^(v)λ^(v)= 0,
    
    		   λ^(v)= 
    		-(Θ^(v))^-1𝐜^(v).

Then λ_a^(v),  𝐠_a^(v) and the next step velocity 𝐯_a(t+δ t)  are solved. 


Using inverse matrix method to solve the Lagrange multipliers is usually considered to be slow. However, no matter how many atoms are in one rotational constraint, the constrained DOF are only 3(6 for roto-translational constraint); we only need to ask for the inverse of two 3×3(6× 6 for roto-translational constraint) matrices. Therefore, using inverse matrix method should be acceptable.

At the next step, we compute all the Lagrange multipliers, constrained positions, and constrained velocities again, then loop until the end of MD. 



§ PROGRAM IMPLEMENTATION AND RESULTS 




 §.§ Program implementation



  §.§.§ Basic usage

We implemented the program in a homebrew version of CP2K v8.2 package<cit.>.

The basic usage of the program can be shown in a sample input file:
[basicstyle=1.,numbers=left, tabsize=2, linewidth= ]
...
   COORD !Cartesian coordinates of system
	N  -0.00000000    0.00000000   -0.09999998  M1
	H   0.00000000    0.95223715    0.23666663  M1
	H   0.82466156   -0.47611858    0.23666663  M1
	H  -0.82466156   -0.47611858    0.23666663  M1
	N   1.61846237    3.12648673    1.12207349  M2
	H   1.61846237    4.07872389    1.45874010  M2
	H   2.44312393    2.65036816    1.45874010  M2
	H   0.79380080    2.65036816    1.45874010  M2
	N  -0.94386640   -1.51381415    2.30701324  M3
	H  -0.94386640   -0.56157700    2.64367985  M3
	H  -0.11920484   -1.98993273    2.64367985  M3
	H  -1.76852797   -1.98993273    2.64367985  M3
   END
...	
   CONSTRAINT !Constraint settings
	CONSTRAINT_INIT T
	   ROTC !Rotational/roto-translational constraint settings
		MOLNAME M1
		ATOMS 1 2 3 4
		POSITIONS0  [ANGSTROM] -0.00000000    0.00000000   -0.09999998  0.00000000    0.95223715    0.23666663  0.82466156   -0.47611858    0.23666663  -0.82466156   -0.47611858    0.23666663
	   END ROTC
   END CONSTRAINT

In this file we fixed the rotation of one ammonia molecule which named as M1, and the other two ammonia molecules M2 and M3 are free of constraints.
Only related text in the input file is shown. Different input blocks are divided by symbol &. Lines 2–15 are coordinates of the system. The first column is the elements of atoms, the followed three columns are the Cartesian coordinates in Angstrom. The last column is the molecule names which are defined as M1, M2, and M3 in the file. Lines 17–24 are settings for the constraints. Lines 19–23 correspond to the “ROTC” block which controls the rotational constraints. “MOLNAME M1” indicates the name which will be constrained is M1. “ATOMS 1 2 3 4” indicates atoms 1–4 in M1 will be constrained. “POSITIONS0 ...” is the rotational reference coordinates at Cartesian coordinates in sequence (x_1, y_1, z_1, x_2, y_2, z_2, x_3, y_3, z_3, ...), the unit is Angstrom. “POSITIONS0 ...” will be transferred to the COM coordinates automatically in the program, which corresponds to the rotational reference coordinates (X_1, Y_1, Z_1, X_2, Y_2, Z_2, …, X_a, Y_a, Z_a) in the Eq.(<ref>).
The usage of  roto-translational constraint algorithm is same as rotational constraint, the algorithm is implemented in another homebrew version of CP2K v8.2.
 




  §.§.§ Center-of-mass definition with periodic boundary condition

Users should note that the definition of the rotational reference center is not unique when using periodic boundary conditions.  As shown in Fig. <ref>, the box in the figure represents the simulation box with the periodic boundary, A and B denote two atoms close to the boundary of the box, and A^' denotes the image of atom A due to periodicity. If the size of the molecule AB is comparable to the size of the box, as in Fig. <ref>(a), its center of mass, i.e., the rotational reference center, will be  at a point O halfway between A and B. If the molecule AB is smaller than half the length of the box, as in Fig. <ref>(b), the rotational reference center will be at the point O^' between A^' and B. Even if we directly define the two atoms A and B as belonging to the same molecule in the input file, we still cannot avoid this ambiguity, so we use Cartesian coordinates to distinguish the two different cases.  In CP2K, as long as “CENTER_COORDINATES” is set as False or default, the coordinates in the input file will be entered as they are, then we can deal with this question automatically by the program. 
If the molecule corresponds to the case in Fig. <ref>(a), the coordinates of atoms A and B can be defined directly, on the contrary, if the molecule corresponds to the case in Fig. <ref>(b), then the coordinates of A^' and B should be defined, instead of the coordinates of A and B.




  §.§.§ Computational details

The running details of the MD trajectory results are as follows: 
The system consists of three ammonia molecules, as in Fig. <ref>.
The HFX (Hartree-Fock exchange) basis set<cit.> and the GTH pseudopotentials<cit.> are used in the AIMD simulation. The simulation box is a cube of 10 Å side length, periodic boundary condition is used. In force field molecular dynamics(FFMD), Amber format force field<cit.> and 25 Å side length box is used for the modeling. 
The initial temperature is 300 K. 
The time step is 1 femtosecond, and the length of each trajectory is 2 picoseconds.
The NVE and the NVT ensemble are used. The NOSÉ thermostator<cit.> is used in the NVT ensemble, and the temperature is set as 300 K. 
 Only AIMD results are shown below. FFMD results are in the Supporting Information. All MD results ran with homebrew version codes of CP2K v8.2.



 §.§ AIMD results



  §.§.§ Rotational constraint results

The snapshot of NVE ensemble AIMD trajectory using rotational constraints are shown in Fig. <ref>(a)(t = 0 fs) and Fig. <ref>(b)(t = 320 fs).
One can check animate graphs in the Supporting Information for more direct view. The result shows the molecule M3 moves normally, and the other two constrained molecules M1 and M2 do not rotate while vibrating and translating, which means our program meets the expectation. We define rotational constraint error as the norm of deviation from Eq. (<ref>)

    σ =  ∑^n_a_a=1𝐫^0_a × m_a 𝐫_a(t) .

In AIMD the max error is σ_max = 9.43 × 10^-7 Å, and in FFMD the max error is σ_max = 1.52 × 10^-6 Å, which are close to numerical errors.


One import question is the numerical stability of the program and insuring the errors will not accumulate during the running of MD. We statistic the fluctuations of conserved quantity, potential and temperature in trajectories, then compare the data with MD without constraints. The results of AIMD are shown in Fig. <ref>–<ref>. 
More results of comparison of FFMD with/without constraints are in Supporting Information. 
The constant quantity(Cons. Qty.) indicates the number stability of the MD.


It corresponds to the total energy in NVE ensemble, and in the NVT ensemble, it corresponds to the total energy of the system containing DOF of the heat bath. The difference of conserved quantity with/without constraints comes from the DOF change by constraints. The results show that no significant difference in the fluctuation size of conserved quantity, potential and temperature.


We also computed the drift and fluctuation of conserved quantity by least-squares fitting. We performed a least-squares fit of data to a straight line respectively. The drift is defined as the difference of the fit at the first and last point. The fluctuation is defined as the root mean square deviation(RMSD) around the least-squares fit. All drift/fluctuation difference from the no constraint MD is showed in Table <ref>. 
Most drift and fluctuation difference values are little than 1×10^-4 a.u. The drift difference in AIMD using NVT ensemble showed a little bigger result with 2.57×10^-4 a.u., which may caused by the collision of two ammonia molecules in the trajectory.




  §.§.§ Roto-translational constraint results

The snapshot of NVE ensemble AIMD trajectory using roto-translational constraints are shown in Fig. <ref>(a)(t = 0 fs) and Fig. <ref>(c)(t = 320 fs).
One can check animate graphs in the Supporting Information for more direct view. The result shows the molecule M3 moves normally, and the other two constrained molecules M1 and M2 do not rotate or translate while vibrating, which means our program meets the expectation. We define the rotational constraint error same as in Eq. (<ref>). In AIMD the max rotational error is σ_max = 3.54 × 10^-7 Å, and in FFMD the max rotational error is σ_max = 4.83 × 10^-6 Å. 
 The translational error is defined as τ = 𝐂^0 - 𝐂(t), which is the vector norm of change of center-of-mass 𝐂. In AIMD the max translational error is τ_max = 1.00 × 10^-6 Å, and in FFMD the max translational error is τ_max = 5.04 × 10^-7 Å, , which are close to numerical error. 
The results of statistic fluctuation is similar with rotational constraint, which is showed in Fig. <ref>–<ref>, no significant difference from the regular MD. The drift and fluctuation difference values are in Table <ref>. The results shows that the roto-translatonal constraint has a good numerical stability.





There is a small difference in initial temperature between constrained MD and regular MD. We set the initial temperature to 300 K in all MD results, when using the roto-translational constraint, the initial temperature is always a little higher than our settings, in Fig. <ref>–<ref> is about 336K. To explain this, we need to understand what is actually happening in the program. Typically, the constraint algorithm does nothing at the zeroth step of MD.
After initializing the atomic velocity according to the set temperature, usually the atomic velocity of the system does not satisfy the constraint condition. Then we apply the constraint to the system, which will do work on the system and cause a change in the total energy. Once the constraint condition is satisfied at a certain step, the constraint algorithm will theoretically do no work on the system at subsequent steps in the MD. This will be shown as a leap in total energy only from zeroth to the first step. 

In order to avoid the total energy leap during the statistic, we used the 
 “CONSTRAINT_INIT T” option, which forces the constraint algorithm to be applied to the zeroth step. Then the program scales the atomic velocities according to the set temperature. 
The overall centre-of-mass velocity of the system is usually zero when scaling, but with the rotational-translational constraint, we fix the translational motion of two molecules, and the remaining molecules will be attracted to or moved away from each other by the interactions, resulting in a total centre-of-mass velocity for the system. Then the CP2K program will also scale the total centre-of-mass velocity, so that the initial temperature will always be different from the temperature we set. We have not modified this problem in order not to affect the use of other functions in CP2K. Although the initial conditions are changed, this problem does not affect the use of the thermostater. Users can also turn off the "CONSTRAINT_INIT" option to avoid problems caused by speed scaling.


§ CONCLUSION



To aid in the study of rotational dynamics in condensed matter systems, we developed a rotational/roto-translational constraint algorithm based on the Velocity Verlet integrator. This algorithm was implemented in a custom version of the CP2K package.

The user interface is friendly. The rotational/roto-translational constraint can be used like other constraint functions in CP2K.

We evaluated the constraint algorithm by AIMD/FFMD using both the NVE and NVT ensembles, and computed the drift and fluctuation of conserved quantity in simulation. The results show that the algorithm can selectively constrain molecules and support stable running of MD.





The algorithm is useful for studying the paddle-wheel mechanism in perovskite materials, phase transitions, and other important rotation-related dynamic problems. In the future, we plan to apply it to the study of solid  electrolytes system, as well as to improve the code structure and the user interface.


	
	This research was sponsored by the National Natural Science Foundation of China (grants 22073035, 22372068, and 21533003), and Graduate Innovation Fund of Jilin University. We also thank the Program for the JLU Science and Technology Innovative Research Team.


	
	
	
	
	
		
		Theoretical derivation of roto-translational constraint, fluctuation of quantities with FFMD, and animates of MD trajectory results. The source code can be found in 
https://github.com/jtyang-chem/rconst_cp2k

		
			
  * SI.pdf: 
			Theoretical derivation of roto-translational constraint.
			Fluctuation of quantities with FFMD.
			
  * animates.zip: Animates of MD trajectory results.
		
		
	
	
	
	
	
	
	
	
	
	
