diff --git a/examples/mdi/README b/examples/mdi/README index 0c4c931519..2dd5f2375e 100644 --- a/examples/mdi/README +++ b/examples/mdi/README @@ -56,13 +56,13 @@ Run with TCP: 3 procs + 4 procs Run with MPI: 1 proc each -% mpirun -np 1 lmp_mpi -mdi "-name driver -role DRIVER -method MPI" -log log.aimd.driver -in in.aimd.driver : -np 1 ../lammps/git/src/lmp_mpi -mdi "-name LAMMPS -role ENGINE -method MPI" -log log.aimd.engine -in in.aimd.engine +% mpirun -np 1 lmp_mpi -mdi "-name driver -role DRIVER -method MPI" -log log.aimd.driver -in in.aimd.driver : -np 1 lmp_mpi -mdi "-name LAMMPS -role ENGINE -method MPI" -log log.aimd.engine -in in.aimd.engine --- Run with MPI: 3 procs + 4 procs -% mpirun -np 3 lmp_mpi -mdi "-name driver -role DRIVER -method MPI" -log log.aimd.driver -in in.aimd.driver : -np 4 ../lammps/git/src/lmp_mpi -mdi "-name LAMMPS -role ENGINE -method MPI" -log log.aimd.engine -in in.aimd.engine +% mpirun -np 3 lmp_mpi -mdi "-name driver -role DRIVER -method MPI" -log log.aimd.driver -in in.aimd.driver : -np 4 lmp_mpi -mdi "-name LAMMPS -role ENGINE -method MPI" -log log.aimd.engine -in in.aimd.engine ------------------------------------------------- ------------------------------------------------- diff --git a/src/MDI/fix_mdi_aimd.cpp b/src/MDI/fix_mdi_aimd.cpp index a998844af7..42b06b31bd 100644 --- a/src/MDI/fix_mdi_aimd.cpp +++ b/src/MDI/fix_mdi_aimd.cpp @@ -185,12 +185,8 @@ void FixMDIAimd::post_force(int vflag) ierr = MDI_Send(buf3all,3*atom->natoms,MDI_DOUBLE,mdicomm); if (ierr) error->all(FLERR,"MDI: >COORDS data"); - // trigger engine to evaluate forces,energy,pressure for current system - - ierr = MDI_Send_command("EVAL",mdicomm); - if (ierr) error->all(FLERR,"MDI: EVAL command"); - // request forces from MDI engine + // this triggers engine to evaluate forces,energy,stress for current system ierr = MDI_Send_command("all(FLERR,"MDI: all(FLERR,"MDI: all(FLERR,"MDI: all(FLERR,"MDI: all(FLERR,"MDI: all(FLERR,"MDI: all(FLERR,"MDI: all(FLERR,"MDI: all(FLERR,"MDI: xprd * domain->yprd * domain->zprd; @@ -283,11 +279,12 @@ void FixMDIAimd::reallocate() void FixMDIAimd::unit_conversions() { - double angstrom_to_bohr,kelvin_to_hartree,ev_to_hartree; + double angstrom_to_bohr,kelvin_to_hartree,ev_to_hartree,second_to_aut; MDI_Conversion_factor("angstrom","bohr",&angstrom_to_bohr); MDI_Conversion_factor("kelvin_energy","hartree",&kelvin_to_hartree); MDI_Conversion_factor("electron_volt","hartree",&ev_to_hartree); + MDI_Conversion_Factor("second","atomic_unit_of_time",&second_to_aut); // length units @@ -312,7 +309,7 @@ void FixMDIAimd::unit_conversions() mdi2lmp_energy = 1.0 / ev_to_hartree; } - // force units + // force units = energy/length mdi2lmp_force = 1.0; lmp2mdi_force = 1.0; @@ -325,8 +322,31 @@ void FixMDIAimd::unit_conversions() mdi2lmp_force = angstrom_to_bohr / ev_to_hartree; } - // pressure units + // pressure or stress units = force/area = energy/volume mdi2lmp_pressure = 1.0; lmp2mdi_pressure = 1.0; + + if (lmpunits == REAL) { + lmp2mdi_pressure = (kelvin_to_hartree / force->boltz) / + (angstrom_to_bohr * angstrom_to_bohr * angstrom_to_bohr) / force->nktv2p; + mdi2lmp_pressure = 1.0 / lmp2mdi_pressure; + } else if (lmpunits == METAL) { + lmp2mdi_pressure = ev_to_hartree / + (angstrom_to_bohr * angstrom_to_bohr * angstrom_to_bohr) / force->nktv2p; + mdi2lmp_pressure = 1.0 / lmp2mdi_pressure; + } + + // velocity units = distance/time + + mdi2lmp_velocity = 1.0; + lmp2mdi_velocity = 1.0; + + if (lmpunits == REAL) { + lmp2mdi_velocity = angstrom_to_bohr / (1.0e-15 * second_to_aut); + mdi2lmp_velocity = 1.0 / lmp2mdi_velocity; + } else if (lmpunits == METAL) { + lmp2mdi_velocity = angstrom_to_bohr / (1.0e-12 * second_to_aut); + mdi2lmp_velocity = 1.0 / lmp2mdi_velocity; + } } diff --git a/src/MDI/fix_mdi_aimd.h b/src/MDI/fix_mdi_aimd.h index 63db7266cf..4cf2b536b0 100644 --- a/src/MDI/fix_mdi_aimd.h +++ b/src/MDI/fix_mdi_aimd.h @@ -54,7 +54,7 @@ class FixMDIAimd : public Fix { double lmp2mdi_energy,mdi2lmp_energy; double lmp2mdi_force,mdi2lmp_force; double lmp2mdi_pressure,mdi2lmp_pressure; - double lmp2mdi_virial,mdi2lmp_virial; + double lmp2mdi_velocity,mdi2lmp_velocity; // buffers for MDI comm diff --git a/src/MDI/mdi_engine.cpp b/src/MDI/mdi_engine.cpp index d089421612..db0636b33f 100644 --- a/src/MDI/mdi_engine.cpp +++ b/src/MDI/mdi_engine.cpp @@ -150,6 +150,7 @@ void MDIEngine::mdi_engine(int narg, char **arg) ibuf1 = ibuf1all = nullptr; maxbuf = 0; + need_evaluation = 1; nbytes = -1; create_atoms_flag = 0; @@ -300,6 +301,7 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) } else if (strcmp(command,">COORDS") == 0) { receive_double3(COORD,0); + need_evaluation = 1; } else if (strcmp(command,">FORCES") == 0) { receive_double3(FORCE,0); @@ -328,11 +330,17 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) send_double3(COORD); } else if (strcmp(command,"ago < 0) { + + update->whichflag = 1; + lmp->init(); + update->integrate->setup(1); + update->whichflag = 0; + + } else { + + // insure potential energy and virial are tallied on this step + + update->ntimestep++; + pe->addstep(update->ntimestep); + press->addstep(update->ntimestep); + + int nflag = neighbor->decide(); + if (nflag == 0) { + comm->forward_comm(); + update->integrate->setup_minimal(0); + modify->clearstep_compute(); + output->thermo->compute(1); + modify->addstep_compute(update->ntimestep+1); + } else { + update->integrate->setup_minimal(1); + modify->clearstep_compute(); + output->thermo->compute(1); + modify->addstep_compute(update->ntimestep+1); + } + } +} + // ---------------------------------------------------------------------- // ---------------------------------------------------------------------- // responses to standard MDI driver commands @@ -1356,57 +1418,6 @@ void MDIEngine::infile() delete [] infile; } -/* ---------------------------------------------------------------------- - EVAL command - compute forces, energy, pressure of current system - can be called multiple times by driver - for a system that is continuously evolving - distinguishes between: - (1) first-time call - (2) system needs reneighboring - (3) system does not need reneighboring - this method does NOT increment timestep ----------------------------------------------------------------------- */ - -void MDIEngine::evaluate() -{ - // NOTE: ago is not a good test - // caller needs to distinguish - // currently cannot call it interspersed with "delete atoms" calls - // b/c whichflag stays set - // separate issue, need to unset whichflag when done - - if (neighbor->ago < 0) { - - update->whichflag = 1; - lmp->init(); - update->integrate->setup(1); - update->whichflag = 0; - - } else { - - // insure potential energy and virial are tallied on this step - - update->ntimestep++; - pe->addstep(update->ntimestep); - press->addstep(update->ntimestep); - - int nflag = neighbor->decide(); - if (nflag == 0) { - comm->forward_comm(); - update->integrate->setup_minimal(0); - modify->clearstep_compute(); - output->thermo->compute(1); - modify->addstep_compute(update->ntimestep+1); - } else { - update->integrate->setup_minimal(1); - modify->clearstep_compute(); - output->thermo->compute(1); - modify->addstep_compute(update->ntimestep+1); - } - } -} - /* ---------------------------------------------------------------------- RESET_BOX command wrapper on library reset_box() method diff --git a/src/MDI/mdi_engine.h b/src/MDI/mdi_engine.h index 36a6fbbf8c..25ed98bb5b 100644 --- a/src/MDI/mdi_engine.h +++ b/src/MDI/mdi_engine.h @@ -55,6 +55,8 @@ class MDIEngine : public Command { class Minimize *minimizer; class Compute *ke,*pe,*press; + int need_evaluation; // 1 if system has changed, else 0 + int nbytes; // NBYTES command value used by other commands // create_atoms state @@ -90,6 +92,8 @@ class MDIEngine : public Command { void mdi_md(); void mdi_optg(); + void evaluate(); + void receive_natoms(); void send_natoms(); void send_ntypes(); @@ -115,7 +119,6 @@ class MDIEngine : public Command { void single_command(); void many_commands(); void infile(); - void evaluate(); void reset_box(); void create_atoms(int); void send_stress();