diff --git a/examples/mdi/series_driver.py b/examples/mdi/series_driver.py index 682a6365ec..6addedcb71 100644 --- a/examples/mdi/series_driver.py +++ b/examples/mdi/series_driver.py @@ -150,9 +150,9 @@ random.seed(seed) for icalc in range(ncalc): - # delete all atoms so can run a new calculation + # define a new system - send_command("delete_atoms group all") + mdi.MDI_Send_command("@INIT_SYS",mdicomm) # define simulation box @@ -166,8 +166,8 @@ for icalc in range(ncalc): # send simulation box to engine - vec = [xlo,ylo,zlo,xhi,yhi,zhi,0.0,0.0,0.0] - mdi.MDI_Send_command("RESET_BOX",mdicomm) + vec = [xhi-xlo,0.0,0.0] + [0.0,yhi-ylo,0.0] + [0.0,0.0,zhi-zlo] + mdi.MDI_Send_command(">CELL",mdicomm) mdi.MDI_Send(vec,9,mdi.MDI_DOUBLE,mdicomm) # create atoms on perfect lattice @@ -203,36 +203,43 @@ for icalc in range(ncalc): # send atoms and their properties to engine - mdi.MDI_Send_command("CREATE_ATOM",mdicomm) + mdi.MDI_Send_command(">NATOMS",mdicomm) mdi.MDI_Send(natoms,1,mdi.MDI_INT,mdicomm) - mdi.MDI_Send_command("CREATE_TYPE",mdicomm) + mdi.MDI_Send_command(">TYPES",mdicomm) mdi.MDI_Send(atypes,natoms,mdi.MDI_INT,mdicomm) - mdi.MDI_Send_command("CREATE_X",mdicomm) + mdi.MDI_Send_command(">COORDS",mdicomm) mdi.MDI_Send(coords,3*natoms,mdi.MDI_DOUBLE,mdicomm) - mdi.MDI_Send_command("CREATE_V",mdicomm) + mdi.MDI_Send_command(">VELOCITIES",mdicomm) mdi.MDI_Send(vels,3*natoms,mdi.MDI_DOUBLE,mdicomm) - mdi.MDI_Send_command("CREATE_GO",mdicomm) # eval or run or minimize + mdi.MDI_Send_command("@EVAL",mdicomm) + if mode == "eval": - mdi.MDI_Send_command("EVAL",mdicomm) + mdi.MDI_Send_command(">STEPS",mdicomm) + mdi.MDI_Send(0,1,mdi.MDI_INT,mdicomm) elif mode == "run": - send_command("run %d" % nsteps) + mdi.MDI_Send_command(">STEPS",mdicomm) + mdi.MDI_Send(steps,1,mdi.MDI_INT,mdicomm) elif mode == "min": - send_command("minimize %g %g 1000 1000" % (tol,tol)) + mdi.MDI_Send_command(">TOLERANCE",mdicomm) + params = [1.0e-4,1.0e-4,1000.0,1000.0] + mdi.MDI_Send(params,4,mdi.MDI_DOUBLE,mdicomm) - # request energy + mdi.MDI_Send_command("@DEFAULT",mdicomm) - mdi.MDI_Send_command("all(FLERR,"MDI: MDI engine is already performing a simulation"); mode = SYS; node_match = false; + } else if (strcmp(command,"@EVAL") == 0) { + if (mode != DEFAULT) + error->all(FLERR,"MDI: MDI engine is already performing a simulation"); + mode = EVAL; + node_match = false; } else if (strcmp(command,"NBYTES") == 0) { nbytes_command(); @@ -466,7 +470,11 @@ int MDIEngine::execute_command(const char *command, MDI_Comm mdicomm) } else if (strcmp(command,"COMMANDS") == 0) { many_commands(); } else if (strcmp(command,"INFILE") == 0) { - infile(); + infile();w + } else if (strcmp(command,">STEPS") == 0) { + send_steps(); + } else if (strcmp(command,">TOLERANCE") == 0) { + send_tolerance(); } else if (strcmp(command,"NATOMS"); MDI_Register_command("@INIT_SYS", ">TYPES"); MDI_Register_command("@INIT_SYS", ">VELOCITIES"); + MDI_Register_command("@INIT_SYS", "@DEFAULT"); MDI_Register_command("@INIT_SYS", "EXIT"); + // custom nodes added by LAMMPS + + MDI_Register_node("@EVAL"); + MDI_Register_command("@EVAL", ">STEPS"); + MDI_Register_command("@EVAL", ">TOLERANC"); + // node for setting up and running a dynamics simulation MDI_Register_node("@INIT_MD"); @@ -749,72 +764,129 @@ void MDIEngine::mdi_optg() } /* ---------------------------------------------------------------------- - initialize a new simulation for single-point calc or dynamics or min + initialize a new simulation ---------------------------------------------------------------------- */ void MDIEngine::mdi_sys() { // engine is now at @INIT_SYS node // receive commands driver sends to define the system - // GO command will force return from INIT_SYS + // another @ command will trigger setup of the system + + sys_natoms_flag = sys_types_flag = sys_charges_flag = 0; + sys_coords_flag = sys_velocities_flag = 0; + sys_cell_flag = sys_cell_displ_flag = 0; engine_node("@INIT_SYS"); - if (strcmp(mdicmd,"@DEFAULT") == 0 || strcmp(mdicmd,"EXIT") == 0) return; + if (strcmp(mdicmd,"EXIT") == 0) return; + + // >CELL, >NATOMS, >TYPES, >COORDS commands must have been issued + + if (sys_cell_flag == 0 || sys_natoms_flag == 0) + error->all(FLERR,"@INIT_SYS requires >NATOMS and >CELL MDI commands"); + if (sys_types_flag == 0 || sys_coords_flag == 0) + error->all(FLERR,"@INIT_SYS requires >TYPES and >COORDS MDI commands"); // clear system via delete_atoms command + + lmp->input->one("delete_atoms group all"); + // lib->reset_box() - // lib->create_atoms() - // then init LAMMPS for new system for MD or min - //lmp->input->command(); + double boxlo[3],boxhi[3]; + double xy,yz,xz; - // initialize a new simulation - - update->whichflag = 1; - timer->init_timeout(); - update->nsteps = 1; - update->ntimestep = 0; - update->firststep = update->ntimestep; - update->laststep = update->ntimestep + update->nsteps; - update->beginstep = update->firststep; - update->endstep = update->laststep; - - lmp->init(); - - // setup the MD simulation - - update->integrate->setup(1); - - // run MD one step at a time until driver sends @DEFAULT or EXIT - // driver can communicate with LAMMPS within each timestep - // by sending a node command which matches a method in FixMDIEngine - - while (true) { - update->whichflag = 1; - timer->init_timeout(); - update->nsteps += 1; - update->laststep += 1; - update->endstep = update->laststep; - output->next = update->ntimestep + 1; - - // single MD timestep - - update->integrate->run(1); - - // driver triggers end of MD loop by sending @DEFAULT or EXIT - - if (strcmp(mdicmd,"@DEFAULT") == 0 || strcmp(mdicmd,"EXIT") == 0) { - modify->delete_fix("MDI_ENGINE_INTERNAL"); - return; - } + if (sys_cell_displ_flag) { + boxlo[0] = sys_cell_displ[0]; + boxlo[1] = sys_cell_displ[1]; + boxlo[2] = sys_cell_displ[2]; + } else { + boxlo[0] = boxlo[1] = boxlo[2] = 0.0; } - // check that return is NOT GO + boxhi[0] = boxlo[0] + sys_cell[0]; + boxhi[1] = boxlo[1] + sys_cell[4]; + boxhi[2] = boxlo[2] + sys_cell[8]; + + xy = sys_cell[3]; + yz = sys_cell[7]; + xz = sys_cell[6]; - engine_node("@INIT_SYS"); + lammps_reset_box(lmp,boxlo,boxhi,xy,yz,xz); - if (strcmp(mdicmd,"@DEFAULT") == 0 || strcmp(mdicmd,"EXIT") == 0) return; + // lib->create_atoms() + // optionally set charges if specified by ">CHARGES" + + if (sys_velocities_flag) + int natoms = lammps_create_atoms(lmp,sys_natoms,NULL,sys_types, + sys_coords,sys_velocities,NULL,1); + else + int natoms = lammps_create_atoms(lmp,sys_natoms,NULL,sys_types, + sys_coords,NULL,NULL,1); + + if (sys_charges_flag) lammps_scatter_atoms(lmp,(char *) "q",1,1,sys_charges); + + // new system + + update->ntimestep = 0; +} + +/* ---------------------------------------------------------------------- + perform a single-point calc or dynamics or min +---------------------------------------------------------------------- */ + +void MDIEngine::mdi_eval() +{ + // engine is now at @EVAL node + // receive commands driver sends to define the evaluation + // another @ command will trigger the calculation to be performed + + engine_node("@EVAL"); + + if (strcmp(mdicmd,"EXIT") == 0) return; + + if (evalmode == POINT or evalmode == MD) { + update->whichflag = 1; + timer->init_timeout(); + + update->nsteps = nsteps; + update->beginstep = update->firststep = update->ntimestep; + update->endstep = update->laststep = update->ntimestep + update->nsteps; + + lmp->init(); + update->integrate->setup(1); + + timer->init(); + timer->barrier_start(); + update->integrate->run(nsteps); + timer->barrier_stop(); + + update->integrate->cleanup(); + } + + if (evalmode == MINIMIZE) { + update->etol = etol[0]; + update->ftol = ftol[1]; + update->nsteps = nsteps; + update->max_eval = max_eval; + + update->whichflag = 2; + timer->init_timeout(); + + update->beginstep = update->firststep = update->ntimestep; + update->endstep = update->laststep = update->firststep + update->nsteps; + + lmp->init(); + update->minimize->setup(); + + timer->init(); + timer->barrier_start(); + update->minimize->run(update->nsteps); + timer->barrier_stop(); + + update->minimize->cleanup(); + } } /* ---------------------------------------------------------------------- @@ -1226,6 +1298,93 @@ void MDIEngine::receive_double3(int which, int addflag) // ---------------------------------------------------------------------- // ---------------------------------------------------------------------- +/* ---------------------------------------------------------------------- + boxhi[0] - domain->boxlo[0]; + celldata[1] = 0.0; + celldata[2] = 0.0; + celldata[3] = domain->xy; + celldata[4] = domain->boxhi[1] - domain->boxlo[1]; + celldata[5] = 0.0; + celldata[6] = domain->xz; + celldata[7] = domain->yz; + celldata[8] = domain->boxhi[2] - domain->boxlo[2]; + + for (int icell = 0; icell < 9; icell++) + celldata[icell] *= lmp2mdi_length; + + int ierr = MDI_Send(celldata,9,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: boxlo[0]; + celldata[1] = domain->boxlo[1]; + celldata[2] = domain->boxlo[2]; + + for (int icell = 0; icell < 3; icell++) + celldata[icell] *= lmp2mdi_length; + + int ierr = MDI_Send(celldata,3,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: compute_scalar(); + double kinetic_energy = ke->compute_scalar(); + double total_energy = potential_energy + kinetic_energy; + total_energy *= lmp2mdi_energy; + + int ierr = MDI_Send(&total_energy,1,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: natoms * MDI_LABEL_LENGTH]; + memset(labels,' ',atom->natoms * MDI_LABEL_LENGTH); + + // NOTE: this loop will not work in parallel + + for (int iatom = 0; iatom < atom->natoms; iatom++) { + std::string label = std::to_string(atom->type[iatom]); + int label_len = std::min(int(label.length()), MDI_LABEL_LENGTH); + strncpy(&labels[iatom * MDI_LABEL_LENGTH], label.c_str(), label_len); + } + + int ierr = MDI_Send(labels,atom->natoms*MDI_LABEL_LENGTH,MDI_CHAR,mdicomm); + if (ierr) error->all(FLERR,"MDI: all(FLERR,"MDI: compute_scalar(); + potential_energy *= lmp2mdi_energy; + int ierr = MDI_Send(&potential_energy,1,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: compute_vector(); + for (int i = 0; i < 6; i++) + vtensor[i] = press->vector[i] * lmp2mdi_pressure; - - - - - - + int ierr = MDI_Send(vtensor,6,MDI_DOUBLE,mdicomm); + if (ierr) error->all(FLERR,"MDI: all(FLERR,"MDI: natoms * MDI_LABEL_LENGTH]; - memset(labels,' ',atom->natoms * MDI_LABEL_LENGTH); - - // NOTE: this loop will not work in parallel - - for (int iatom = 0; iatom < atom->natoms; iatom++) { - std::string label = std::to_string(atom->type[iatom]); - int label_len = std::min(int(label.length()), MDI_LABEL_LENGTH); - strncpy(&labels[iatom * MDI_LABEL_LENGTH], label.c_str(), label_len); - } - - int ierr = MDI_Send(labels,atom->natoms*MDI_LABEL_LENGTH,MDI_CHAR,mdicomm); - if (ierr) error->all(FLERR,"MDI: compute_scalar(); - double kinetic_energy = ke->compute_scalar(); - double total_energy = potential_energy + kinetic_energy; - total_energy *= lmp2mdi_energy; - - int ierr = MDI_Send(&total_energy,1,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR,"MDI: compute_scalar(); - potential_energy *= lmp2mdi_energy; - - int ierr = MDI_Send(&potential_energy,1,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR,"MDI: compute_scalar(); - kinetic_energy *= lmp2mdi_energy; - - int ierr = MDI_Send(&kinetic_energy,1,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR,"MDI: boxhi[0] - domain->boxlo[0]; - celldata[1] = 0.0; - celldata[2] = 0.0; - celldata[3] = domain->xy; - celldata[4] = domain->boxhi[1] - domain->boxlo[1]; - celldata[5] = 0.0; - celldata[6] = domain->xz; - celldata[7] = domain->yz; - celldata[8] = domain->boxhi[2] - domain->boxlo[2]; - - for (int icell = 0; icell < 9; icell++) - celldata[icell] *= lmp2mdi_length; - - int ierr = MDI_Send(celldata,9,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR,"MDI: boxlo[0]; - celldata[1] = domain->boxlo[1]; - celldata[2] = domain->boxlo[2]; - - for (int icell = 0; icell < 3; icell++) - celldata[icell] *= lmp2mdi_length; - - int ierr = MDI_Send(celldata,3,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR,"MDI: compute_scalar(); + kinetic_energy *= lmp2mdi_energy; - if (atom->natoms > 0) - error->all(FLERR,"MDI RESET_BOX cannot be used when atoms exist"); - - ierr = MDI_Recv(values,9,MDI_DOUBLE,mdicomm); - MPI_Bcast(values,9,MPI_DOUBLE,0,world); - - lammps_reset_box(lmp,&values[0],&values[3],values[6],values[7],values[8]); -} - -/* ---------------------------------------------------------------------- - CREATE_ATOM command - wrapper on library create_atoms() method - requires simulation box be defined - allows caller to define a new set of atoms - with their IDs, types, coords, velocities, image flags - called in stages via flag - since MDI plugin mode only allows 1 MDI Send/Recv per MDI command - assumes current atom->natoms set by >NATOMS command is correct ----------------------------------------------------------------------- */ - -void MDIEngine::create_atoms(int flag) -{ - int ierr; - - // NOTE: error check on imageint = INT - - if (flag == CREATE_ATOM) { - - if (create_atoms_flag) - error->all(FLERR,"MDI CREATE_ATOM already in progress"); - - ierr = MDI_Recv(&create_natoms,1,MDI_INT,mdicomm); - MPI_Bcast(&create_natoms,1,MPI_INT,0,world); - - create_atoms_flag = 1; - create_id = nullptr; - create_type = nullptr; - create_x = nullptr; - create_v = nullptr; - create_image = nullptr; - - } else if (flag == CREATE_ID) { - - if (!create_atoms_flag) error->all(FLERR,"MDI CREATE_ATOM not in progress"); - if (create_id) error->all(FLERR,"MDI CREATE_ATOM already in progress"); - - int natom = create_natoms; - memory->create(create_id,natom,"mdi:create_id"); - ierr = MDI_Recv(create_id,natom,MDI_INT,mdicomm); - MPI_Bcast(create_id,natom,MPI_INT,0,world); - - } else if (flag == CREATE_TYPE) { - - if (!create_atoms_flag) error->all(FLERR,"MDI CREATE_ATOM not in progress"); - if (create_type) error->all(FLERR,"MDI CREATE_ATOM already in progress"); - - int natom = create_natoms; - if (create_type) error->all(FLERR,"MDI CREATE_ATOM already in progress"); - memory->create(create_type,natom,"mdi:create_type"); - ierr = MDI_Recv(create_type,natom,MDI_INT,mdicomm); - MPI_Bcast(create_type,natom,MPI_INT,0,world); - - } else if (flag == CREATE_X) { - - if (!create_atoms_flag) error->all(FLERR,"MDI CREATE_ATOM not in progress"); - if (create_x) error->all(FLERR,"MDI CREATE_ATOM already in progress"); - - int natom = create_natoms; - if (create_x) error->all(FLERR,"MDI CREATE_ATOM already in progress"); - memory->create(create_x,3*natom,"mdi:create_x"); - ierr = MDI_Recv(create_x,3*natom,MDI_DOUBLE,mdicomm); - MPI_Bcast(create_x,3*natom,MPI_DOUBLE,0,world); - - } else if (flag == CREATE_V) { - - if (!create_atoms_flag) error->all(FLERR,"MDI CREATE_ATOM not in progress"); - if (create_v) error->all(FLERR,"MDI CREATE_ATOM already in progress"); - - int natom = create_natoms; - if (create_v) error->all(FLERR,"MDI CREATE_ATOM already in progress"); - memory->create(create_v,3*natom,"mdi:create_x"); - ierr = MDI_Recv(create_v,3*natom,MDI_DOUBLE,mdicomm); - MPI_Bcast(create_v,3*natom,MPI_DOUBLE,0,world); - - } else if (flag == CREATE_IMAGE) { - - if (!create_atoms_flag) error->all(FLERR,"MDI CREATE_ATOM not in progress"); - if (create_image) error->all(FLERR,"MDI CREATE_ATOM already in progress"); - - int natom = create_natoms; - if (create_image) error->all(FLERR,"MDI CREATE_ATOM already in progress"); - memory->create(create_image,natom,"mdi:create_image"); - ierr = MDI_Recv(create_image,natom,MDI_INT,mdicomm); - MPI_Bcast(create_image,natom,MPI_INT,0,world); - - } else if (flag == CREATE_GO) { - - if (!create_atoms_flag) error->all(FLERR,"MDI CREATE_ATOM not in progress"); - if (!create_type || !create_x) - error->all(FLERR,"MDI: CREATE_ATOM requires types and coords"); - - int ncreate = lammps_create_atoms(lmp,create_natoms,create_id,create_type, - create_x,create_v,create_image,1); - - if (ncreate != create_natoms) - error->all(FLERR, "MDI: CREATE ATOM created atoms != sent atoms"); - - // clean up create_atoms state - - create_atoms_flag = 0; - memory->destroy(create_id); - memory->destroy(create_type); - memory->destroy(create_x); - memory->destroy(create_v); - memory->destroy(create_image); - } -} - -/* ---------------------------------------------------------------------- - compute_vector(); - for (int i = 0; i < 6; i++) - vtensor[i] = press->vector[i] * lmp2mdi_pressure; - - int ierr = MDI_Send(vtensor,6,MDI_DOUBLE,mdicomm); - if (ierr) error->all(FLERR,"MDI: all(FLERR,"MDI: