diff --git a/src/Control.h b/src/Control.h index b3401c1d..bd4a7d96 100644 --- a/src/Control.h +++ b/src/Control.h @@ -371,6 +371,8 @@ class Control // 10 or larger means CG, otherwise MG V-cycles bool MGPoissonSolver() { return (diel_flag_ / 10 == 0); } + bool LangevinThermostat() { return (thermostat_type == 1); } + // // data // diff --git a/src/HDFrestart.cc b/src/HDFrestart.cc index 0322bfea..0fac518e 100644 --- a/src/HDFrestart.cc +++ b/src/HDFrestart.cc @@ -457,7 +457,8 @@ HDFrestart::HDFrestart(const std::string& filename, const pb::PEenv& pes, verbosity_ = 0; closed_ = false; - //(*MPIdata::sout)<<"HDFrestart::HDFrestart(), filename="<& data) // send data to inactive PEs if (gather_data_x_) gatherDataXdir(data); + if (useHdf5p()) + { + data.erase(std::remove(data.begin(), data.end(), -1), data.end()); + } + return 0; } @@ -1928,7 +1934,7 @@ int HDFrestart::readAtomicData( std::string datasetname, std::vector& data) { if (onpe0) - (*MPIdata::sout) << "Read ionic positions from hdf5 file" << std::endl; + (*MPIdata::sout) << "Read atomic data from hdf5 file" << std::endl; if (active_) { @@ -1968,6 +1974,10 @@ int HDFrestart::readAtomicData( } } + if (useHdf5p()) + { + data.erase(std::remove(data.begin(), data.end(), 1e+32), data.end()); + } if (gather_data_x_) gatherDataXdir(data); return 0; @@ -2094,6 +2104,11 @@ int HDFrestart::readAtomicData( data.push_back(t); } + if (useHdf5p()) + { + data.erase(std::remove(data.begin(), data.end(), ""), data.end()); + } + return 0; } @@ -2128,7 +2143,6 @@ int HDFrestart::readRestartRandomStates(std::vector& data) dim = (int)H5Dget_storage_size(dataset_id) / sizeof(unsigned short); } - if (dim > 0) { data.resize(dim); @@ -2148,12 +2162,10 @@ int HDFrestart::readRestartRandomStates(std::vector& data) } if (!data.empty()) - if (data[0] != data[0]) + if (std::isnan(data[0])) { MGMOL_HDFRESTART_FAIL( - "ERROR: HDFrestart::readRestartRandomStates() " - "--- data[0]=" - << data[0]); + "readRestartRandomStates() is NaN"); return -2; } } diff --git a/src/Ions.cc b/src/Ions.cc index 913a0f55..e93d82ea 100644 --- a/src/Ions.cc +++ b/src/Ions.cc @@ -911,8 +911,10 @@ void Ions::initFromRestartFile(HDFrestart& h5_file) assert(at_numbers.size() == at_nlprojIds.size()); num_ions_ = at_names.size(); - mmpi.allreduce(&num_ions_, 1, MPI_SUM); - + if (!h5_file.useHdf5p()) + { + mmpi.allreduce(&num_ions_, 1, MPI_SUM); + } if (onpe0 && ct.verbose > 0) { (*MPIdata::sout) << "Ions::setFromRestartFile(), read " << num_ions_ @@ -947,9 +949,21 @@ void Ions::initFromRestartFile(HDFrestart& h5_file) } readRestartPositions(h5_file); readRestartVelocities(h5_file); - readRestartRandomStates(h5_file); + if (ct.LangevinThermostat()) readRestartRandomStates(h5_file); readLockedAtomNames(h5_file); + // remove atoms from local list if not local + for (std::vector::iterator it = local_ions_.begin(); + it != local_ions_.end();) + { + double p[3]; + (*it)->getPosition(p); + if (!inLocalIons(p[0], p[1], p[2])) + it = local_ions_.erase(it); + else + ++it; + } + // rescale all velocities by factor specified in input rescaleVelocities(ct.VelocityScalingFactor()); diff --git a/src/restart.cc b/src/restart.cc index 2fa76ecb..247cca06 100644 --- a/src/restart.cc +++ b/src/restart.cc @@ -117,7 +117,7 @@ int MGmol::write_hdf5(HDFrestart& h5f_file, ions.writeAtomicIDs(h5f_file); ions.writeAtomicNLprojIDs(h5f_file); ions.writePositions(h5f_file); - ions.writeRandomStates(h5f_file); + if (ct.LangevinThermostat()) ions.writeRandomStates(h5f_file); ions.writeVelocities(h5f_file); ions.writeForces(h5f_file); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b29c1aad..baf1c451 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -504,6 +504,14 @@ add_test(NAME testMD_D72 ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/coords.in ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/lrs.in ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +add_test(NAME testHDF5single + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt + ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/md.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/h2o.xyz + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) add_test(NAME testMD_MVP COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MD_MVP/test.py ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} diff --git a/tests/DMandEnergyAndForces/test.py b/tests/DMandEnergyAndForces/test.py index c36f0767..b21e6e72 100755 --- a/tests/DMandEnergyAndForces/test.py +++ b/tests/DMandEnergyAndForces/test.py @@ -3,6 +3,7 @@ import os import subprocess import string +import shutil print("Test DMandEnergyAndForces...") @@ -33,6 +34,8 @@ output = subprocess.check_output(command,shell=True) lines=output.split(b'\n') +shutil.rmtree('WF') + #analyse output energies=[] for line in lines: diff --git a/tests/HDF5single/h2o.xyz b/tests/HDF5single/h2o.xyz new file mode 100644 index 00000000..d5171c8b --- /dev/null +++ b/tests/HDF5single/h2o.xyz @@ -0,0 +1,6 @@ +3 + +O 0.00 0.00 0.00 +H -0.76 0.59 0.00 +H 0.76 0.59 0.00 + diff --git a/tests/HDF5single/md.cfg b/tests/HDF5single/md.cfg new file mode 100644 index 00000000..1ff2adab --- /dev/null +++ b/tests/HDF5single/md.cfg @@ -0,0 +1,34 @@ +verbosity=3 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=48 +ny=48 +nz=48 +[Domain] +ox=-4.5 +oy=-4.5 +oz=-4.5 +lx=9. +ly=9. +lz=9. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=MD +[MD] +num_steps=5 +dt=40. +[Quench] +max_steps=24 +atol=1.e-8 +[Restart] +input_level=4 +input_filename=WF +input_type=single_file +output_level=4 +output_filename=WF_MD +output_type=single_file +[Coloring] +scope=global diff --git a/tests/HDF5single/mgmol.cfg b/tests/HDF5single/mgmol.cfg new file mode 100644 index 00000000..4dba942a --- /dev/null +++ b/tests/HDF5single/mgmol.cfg @@ -0,0 +1,31 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=48 +ny=48 +nz=48 +[Domain] +ox=-4.5 +oy=-4.5 +oz=-4.5 +lx=9. +ly=9. +lz=9. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=QUENCH +[Quench] +max_steps=120 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=1.5 +[Restart] +output_level=4 +output_filename=WF +output_type=single_file +[Coloring] +scope=global diff --git a/tests/HDF5single/test.py b/tests/HDF5single/test.py new file mode 100755 index 00000000..080ee0ba --- /dev/null +++ b/tests/HDF5single/test.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python +import sys +import os +import subprocess +import string + +print("Test test_rho_restart...") + +nargs=len(sys.argv) + +mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3] +for i in range(4,nargs-7): + mpicmd = mpicmd + " "+sys.argv[i] +print("MPI run command: {}".format(mpicmd)) + +mgmol_exe = sys.argv[nargs-5] +input1 = sys.argv[nargs-4] +input2 = sys.argv[nargs-3] +coords = sys.argv[nargs-2] +print("coordinates file: %s"%coords) + +#create links to potentials files +dst1 = 'pseudo.H_ONCV_PBE_SG15' +src1 = sys.argv[-1] + '/' + dst1 + +dst2 = 'pseudo.O_ONCV_PBE_SG15' +src2 = sys.argv[-1] + '/' + dst2 + +if not os.path.exists(dst1): + print("Create link to %s"%dst1) + os.symlink(src1, dst1) + +if not os.path.exists(dst2): + print("Create link to %s"%dst2) + os.symlink(src2, dst2) + +#run mgmol to generate initial ground state +command = "{} {} -c {} -i {}".format(mpicmd,mgmol_exe,input1,coords) +print("Run command: {}".format(command)) + +output = subprocess.check_output(command,shell=True) +lines=output.split(b'\n') + +#run MD +command = "{} {} -c {} -i {}".format(mpicmd,mgmol_exe,input2,coords) +print("Run command: {}".format(command)) +output = subprocess.check_output(command,shell=True) +lines=output.split(b'\n') + +os.remove('WF') + +print("Check energy conservation...") +tol = 1.e-4 +energy = 0. +count = 0 +for line in lines: + if line.count(b'Total') and line.count(b'Energy'): + print(line) + count=count+1 + words=line.split() + + energy=eval(words[2]) + if count==1: + first_energy=energy + + if count>1 and abs(energy-first_energy)>tol: + print("ERROR Energy = {} != {}".format(energy,first_energy)) + sys.exit(1) + +if count<4: + print("ERROR needs 4 energy values for checking conservation!") + sys.exit(1) + +print("Test SUCCESSFUL!") +sys.exit(0) diff --git a/tests/RestartEnergyAndForces/test.py b/tests/RestartEnergyAndForces/test.py index 349434e8..b62d39f8 100755 --- a/tests/RestartEnergyAndForces/test.py +++ b/tests/RestartEnergyAndForces/test.py @@ -3,6 +3,7 @@ import os import subprocess import string +import shutil print("Test RestartEnergyAndForces...") @@ -54,14 +55,14 @@ ref_energy=energy break -#sys.exit(0) - #run test command = "{} {} -c {} -i {}".format(mpicmd,test_exe,input2,coords) print("Run command: {}".format(command)) output = subprocess.check_output(command,shell=True) lines=output.split(b'\n') +shutil.rmtree('WF') + test_energy=1.e18 for line in lines: if line.count(b'%%'): diff --git a/tests/WFEnergyAndForces/test.py b/tests/WFEnergyAndForces/test.py index 45420ddf..6cac8c93 100755 --- a/tests/WFEnergyAndForces/test.py +++ b/tests/WFEnergyAndForces/test.py @@ -3,6 +3,7 @@ import os import subprocess import string +import shutil print("Test WFEnergyAndForces...") @@ -37,6 +38,8 @@ output = subprocess.check_output(command,shell=True) lines=output.split(b'\n') +shutil.rmtree('WF') + #analyse output energies=[] for line in lines: