input files for molecular dynamics calculation using cp.x

+ fix latex eqn. in problem statement page
+ update docusaurus to v3.1.0

docs/theory/problem-statement.md

@@ -54,11 +54,12 @@ However, the challenge is to solve the Schrödinger equation as a real physical
 system is consists of a large number of atoms. The Schrödinger equation becomes
 coupled many-body equation.
 
-$$\left[-\frac{\hbar}{2m} \sum_{i=1}^N \nabla_i^2 + \sum_{i=1}^NV(\textbf{r}_i)
+$$
+\left[-\frac{\hbar}{2m} \sum_{i=1}^N \nabla_i^2 + \sum_{i=1}^NV(\textbf{r}_i)
 + \sum_{i=1}^N \sum_{j<i}U(\textbf{r}_i, \textbf{r}_j)\right]\psi(\textbf{r}_1,
-\textbf{r}_2, ..., \textbf{r}_N)$$
-$$= E\psi(\textbf{r}_1, \textbf{r}_2, ...,
-\textbf{r}_N)$$
+\textbf{r}_2, ..., \textbf{r}_N) = E\psi(\textbf{r}_1, \textbf{r}_2, ...,
+\textbf{r}_N)
+$$
 
 With today's available computing power, it is far from feasible to solve the
 actual electronic wavefunction of a condensed matter system, where $N$ is of the

package.json

@@ -15,8 +15,8 @@
     "write-heading-ids": "docusaurus write-heading-ids"
   },
   "dependencies": {
-    "@docusaurus/core": "3.0.1",
-    "@docusaurus/preset-classic": "3.0.1",
+    "@docusaurus/core": "3.1.0",
+    "@docusaurus/preset-classic": "3.1.0",
     "@easyops-cn/docusaurus-search-local": "^0.40.1",
     "@mdx-js/react": "^3.0.0",
     "@svgr/webpack": "^5.5.0",

src/H2O-md-cp/README.md

@@ -0,0 +1,33 @@
+# Molecular Dynamics (using cp.x)
+
+- First run relaxation calculation. We can use either `pw.x` or `cp.x` for this
+- Add random displacement.
+
+> [!WARNING]
+> `cp.x` exited while using PAW pseudopotential (not implemented), later I switched to USPP (GBRV) pseudopotentials.
+
+```
+mpirun -np 4 cp.x -i cp.h2o_rand_disp.in | tee cp.h2o_rand_disp.out
+```
+
+- For MD calculation, calculate force and stress for DeePMD. Run md calculation.
+```console
+mpirun -np 4 cp.x -i cp.md.h2o.in | tee cp.md.h2o.out
+```
+
+- Post-processing:
+
+```
+gfortran cp2xsf.f90 -o cp2xsf.x
+
+# execute from the outdir of previous calculation
+cp2xsf.x < cp2xsf.in > cp2xsf.out
+```
+
+The output `.asxf` file can be loaded in Xcrysden for visualization.
+
+Alternatively:
+```
+cppp.x < cppp.in > cppp.out
+```
+

src/H2O-md-cp/cp.h2o_rand_disp.in

@@ -0,0 +1,52 @@
+&control
+    calculation   = 'cp'
+    restart_mode  = 'from_scratch'
+    prefix        = 'h2o'
+    outdir        = './tmp/'
+    pseudo_dir    = '../pseudos/'
+    nstep         = 3000
+    iprint        = 10
+    isave         = 100
+    dt            = 5.0,
+    ndr           = 50
+    ndw           = 50
+    tstress       = .true.
+    tprnfor       = .true.
+    etot_conv_thr = 1.d-9
+    ekin_conv_thr = 1.d-7
+/
+
+&system
+    ibrav     = 1
+    celldm(1) = 15
+    nat       = 3
+    ntyp      = 2
+    ecutwfc   = 50
+    ecutrho   = 400
+    nr1b      = 20
+    nr2b      = 20
+    nr3b      = 20
+/
+
+&electrons
+    electron_dynamics = 'damp'
+    electron_damping  = 0.05,
+    emass             = 300.0
+    orthogonalization = 'ortho'
+    ortho_eps         = 1d-11
+/
+
+&ions
+    ion_dynamics = 'none'
+    tranp        = .true.       ! randomize atomic positions
+    amprp        = 0.05         ! displacement amplitude for the randomization
+/
+
+ATOMIC_SPECIES
+ H  1.00784  H_pbe_v1.4.uspp.F.upf
+ O  15.999   O_pbe_v1.2.uspp.F.upf
+
+ATOMIC_POSITIONS (angstrom)
+O        5.000425817   4.781574242   4.999840589
+H        4.999764264   5.408907510   5.779839929
+H        4.999809919   5.409518248   4.220319481

src/H2O-md-cp/cp.md.h2o.in

@@ -0,0 +1,49 @@
+&control
+    calculation  = 'cp'
+    restart_mode = 'reset_counters'
+    prefix       = 'h2o'
+    outdir       = './tmp/'
+    pseudo_dir   = '../pseudos/'
+    nstep        = 1000
+    iprint       = 10
+    isave        = 100
+    dt           = 5.0
+    ndr          = 50
+    ndw          = 51
+    tstress      = .true.
+    tprnfor      = .true.
+/
+
+&system
+    ibrav     = 1
+    celldm(1) = 15
+    nat       = 3
+    ntyp      = 2
+    ecutwfc   = 50
+    ecutrho   = 400
+    nr1b      = 20
+    nr2b      = 20
+    nr3b      = 20
+/
+
+&electrons
+    electron_dynamics   = 'verlet'
+    electron_velocities = 'zero'
+    emass               = 300.0
+    orthogonalization   = 'ortho'
+    ortho_eps           = 1d-11
+/
+
+&ions
+    ion_dynamics = 'verlet'
+    ion_velocities = 'zero'
+/
+
+ATOMIC_SPECIES
+ H  1.00784  H_pbe_v1.4.uspp.F.upf
+ O  15.999   O_pbe_v1.2.uspp.F.upf
+
+ATOMIC_POSITIONS (angstrom)
+O        5.000425817   4.781574242   4.999840589
+H        4.999764264   5.408907510   5.779839929
+H        4.999809919   5.409518248   4.220319481