Restarting QMD calculations gives erroneous trajectories

[dblinger]

Describe the bug
QMD (Gaussian basis set ab initio MD) gives incorrect trajectory upon restarting. Atomic positions and momenta are different for a calculation with, e.g., 50 time steps versus a calculation with 25 steps followed by restart-type calculation of an additional 25 steps.

Describe settings used
Report all the environment variables used:
export NWCHEM_MODULES=smallqm
export NWCHEM_TOP=/home/d54/nwchem-7.0.2
export NWCHEM_TARGET=LINUX64
export USE_MPI=y
export USE_MPIF=y
export USE_MPIF4=y
export ARMCI_NETWORK=OPENIB
#export BUILD_OPENBLAS=y
#export BUILD_SCALAPACK=y
export USE_SCALAPACK=y
export USE_64TO32=n
#MKL scalapack requires static linking as follows:
export SCALAPACK="-L$MKLROOT/lib/intel64 -Wl,-Bstatic -Wl,--start-group -lmkl_scalapack_ilp64 -lmkl_intel_ilp64 -lmkl_core -lmkl_sequential -lmkl_blacs_openmpi_ilp64 -Wl,--end-group -Wl,-Bdynamic -lpthread -lm"
export SCALAPACK_LIB="$SCALAPACK"
export SCALAPACK_SIZE=8
export LAPACK_SIZE=8
export BLAS_SIZE=8
export BLASOPT="-L$MKLROOT/lib/intel64 -lmkl_intel_ilp64 -lmkl_core -lmkl_sequential -lmkl_core -liomp5 -lpthread -lm"
export LAPACK_LIB="-L$MKLROOT/lib/intel64 -lmkl_intel_ilp64 -lmkl_core -lmkl_sequential -lmkl_core -liomp5 -lpthread -lm"

Report what operating system and distribution you are using:
CentOS Linux release 7.9.2009 (Core)

Attach log files
Attach as many log files as possible.

• stdout/stderr of the NWChem execution
  (N/A)

• complete makefile log
  make.log

• $NWCHEM_TOP/src/tools/build/config.log
  config.log

• $NWCHEM_TOP/src/tools/build/comex/config.log
  comex_config.log

• debugging stack
  (N/A)

To Reproduce

1. Steps to reproduce the behavior:

Extract attached tarball.

In the extracted directory, enter the "full_trajectory" subdirectory and run "input.nw" to completion. The trajectory data produced by this calculation (test.xyz) will serve as the ground truth.

Run input file "input.nw" in subdirectory "first_half" to completion. Copy test.* from the "first_half" subdirectory into "second_half". Enter the directory "second_half" and run "input.nw" to completion. This is a restart-type job that will append trajectory information to "test.xyz" from the calculation carried out in "first_half".

Observe the discrepancy between the positions and momenta from the un-restarted calculation in "full_trajectory/test.xyz" and those in "second_half/test.xyz" that result from the restart.

3. Attach all the input files required to run:
   bug_report.tar.gz

Expected behavior
The nuclear positions and momenta should be identical for the restarted and un-restarted QMD calculations.

[edoapra]

@dblinger thank you very much for the detailed bug report.

[edoapra]

@dblinger By any chance, did you try to verify if the bug exists in the current master branch>

[dblinger]

Hi Edo. Sorry, I didn't think to do this before posting. I did just now clone and build the current master branch (commit 9fdbac6) and I see the same behavior as with the 7.0.2 release.

[edoapra]

Thank you very much for doing this additional testing.

[dblinger]

I should probably acknowledge that I understand that for adiabatic QMD calculations, one can always take the positions and momenta from the last timestep of a prior calculation and prepare input files for a subsequent calculation manually to pick up the propagation where it left off. In this regard, the formal restart capability for QMD would seem to be only a convenience for users. However, in the case of NAMD-type calculations there is no path to read in the initial density matrix if not invoking a restart, so such a workaround is unavailable in that case in my understanding.