From 15999f6518249d6b3104c502b1928dda0fca476f Mon Sep 17 00:00:00 2001 From: ljwoods2 Date: Tue, 12 Nov 2024 14:30:10 -0700 Subject: [PATCH 01/37] initial commit --- src/MISC/fix_imd.cpp | 779 +++++++++++++++++++++++++++++++++++++------ src/MISC/fix_imd.h | 26 +- 2 files changed, 704 insertions(+), 101 deletions(-) diff --git a/src/MISC/fix_imd.cpp b/src/MISC/fix_imd.cpp index 41ce97ee65..c295fdfaf8 100644 --- a/src/MISC/fix_imd.cpp +++ b/src/MISC/fix_imd.cpp @@ -76,6 +76,7 @@ negotiate an appropriate license for such distribution." #endif #include +#include using namespace LAMMPS_NS; using namespace FixConst; @@ -358,7 +359,6 @@ typedef struct { } IMDheader; #define IMDHEADERSIZE 8 -#define IMDVERSION 2 typedef enum IMDType_t { IMD_DISCONNECT, /**< close IMD connection, leaving sim running */ @@ -370,7 +370,15 @@ typedef enum IMDType_t { IMD_MDCOMM, /**< MDComm style force data */ IMD_PAUSE, /**< pause the running simulation */ IMD_TRATE, /**< set IMD update transmission rate */ - IMD_IOERROR /**< indicate an I/O error */ + IMD_IOERROR, /**< indicate an I/O error */ + /* IMDv3 only headers */ + IMD_SESSIONINFO, + IMD_RESUME, + IMD_TIME, + IMD_BOX, + IMD_VELOCITIES, + IMD_FORCES, + IMD_WAIT, } IMDType; /**< IMD command message type enumerations */ typedef struct { @@ -386,8 +394,20 @@ typedef struct { float Eimpr; /**< Improper energy, Kcal/mol */ } IMDEnergies; /**< IMD simulation energy report structure */ +/* IMDv3 only */ +typedef struct IMDSessionInfo { + bool time; + bool box; + bool coords; + bool unwrap; + bool velocities; + bool forces; + bool energies; +} IMDSessionInfo; + /** Send control messages - these consist of a header with no subsequent data */ -static int imd_handshake(void *); /**< check endianness, version compat */ +static int imd_handshake_v2(void *); /**< check endianness, version compat */ +static int imd_handshake_v3(void *, IMDSessionInfo *); /** Receive header and data */ static IMDType imd_recv_header(void *, int32 *); /** Receive MDComm-style forces, units are Kcal/mol/angstrom */ @@ -426,19 +446,22 @@ static void imdsock_destroy(void *); * The implementation follows at the end of the file. * ***************************************************************/ -/* struct for packed data communication of coordinates and forces. */ +/* struct for packed data communication of coordinates, velocities, and forces. */ struct commdata { tagint tag; float x,y,z; }; +MPI_Datatype MPI_CommData; + /*************************************************************** * create class and parse arguments in LAMMPS script. Syntax: - * fix ID group-ID imd [unwrap (on|off)] [fscale ] + * fix ID group-ID imd [version (2|3)] [unwrap (on|off)] [fscale ] [time (on|off)] [box (on|off)] [coordinates (on|off)] [velocities (on|off)] [forces (on|off)] ***************************************************************/ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { + if (narg < 4) error->all(FLERR,"Illegal fix imd command"); @@ -447,12 +470,21 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal fix imd parameter: port < 1024"); /* default values for optional flags */ + imd_version = 2; + unwrap_flag = 0; nowait_flag = 0; connect_msg = 1; imd_fscale = 1.0; imd_trate = 1; + /* IMDv3-only flags. Aren't stored as class attributes since they're converted into IMDSessionInfo */ + int time_flag = 1; + int box_flag = 1; + int coord_flag = 1; + int vel_flag = 1; + int force_flag = 1; + /* parse optional arguments */ int iarg = 4; while (iarg+1 < narg) { @@ -464,9 +496,19 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : imd_fscale = utils::numeric(FLERR,arg[iarg+1],false,lmp); } else if (0 == strcmp(arg[iarg], "trate")) { imd_trate = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - } else { - error->all(FLERR,"Unknown fix imd parameter"); - } + } else if (0 == strcmp(arg[iarg], "version")) { + imd_version = utils::inumeric(FLERR,arg[iarg+1],false,lmp); + } else if (0 == strcmp(arg[iarg], "time")) { + time_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); + } else if (0 == strcmp(arg[iarg], "box")) { + box_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); + } else if (0 == strcmp(arg[iarg], "coordinates")) { + coord_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); + } else if (0 == strcmp(arg[iarg], "velocities")) { + vel_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); + } else if (0 == strcmp(arg[iarg], "forces")) { + force_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); + } else error->all(FLERR,"Unknown fix imd parameter"); ++iarg; ++iarg; } @@ -474,12 +516,43 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : if (imd_trate < 1) error->all(FLERR,"Illegal fix imd parameter. trate < 1."); + if (imd_version != 2 && imd_version != 3) + error->all(FLERR,"Illegal fix imd parameter. version != 2 or 3."); + + imdsinfo = new IMDSessionInfo; + + /* In IMDv2 in LAMMPS, only coordinates are sent*/ + if (imd_version == 2) { + imdsinfo->time = false; + imdsinfo->box = false; + imdsinfo->coords = true; + imdsinfo->unwrap = unwrap_flag; + imdsinfo->velocities = false; + imdsinfo->forces = false; + imdsinfo->energies = false; + } + + if (imd_version == 3) { + imdsinfo->time = time_flag; + imdsinfo->box = box_flag; + imdsinfo->coords = coord_flag; + imdsinfo->unwrap = unwrap_flag; + imdsinfo->velocities = vel_flag; + imdsinfo->forces = force_flag; + imdsinfo->energies = false; + } + + bigint n = group->count(igroup); if (n > MAXSMALLINT) error->all(FLERR,"Too many atoms for fix imd"); num_coords = static_cast (n); - MPI_Comm_rank(world,&me); + // Define MPI dtype for passing x/v/f data + MPI_Type_contiguous(4, MPI_FLOAT, &MPI_CommData); + MPI_Type_commit(&MPI_CommData); + MPI_Comm_rank(world,&me); + /* initialize various imd state variables. */ clientsock = nullptr; localsock = nullptr; @@ -487,14 +560,33 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : imd_inactive = 0; imd_terminate = 0; imd_forces = 0; - force_buf = nullptr; + recv_force_buf = nullptr; maxbuf = 0; - msgdata = nullptr; - msglen = 0; - comm_buf = nullptr; + coord_data = nullptr; + vel_data = nullptr; + force_data = nullptr; idmap = nullptr; rev_idmap = nullptr; + + msglen = 0; + if (imdsinfo->time) { + msglen += 24+IMDHEADERSIZE; + } + if (imdsinfo->box) { + msglen += 9*4+IMDHEADERSIZE; + } + if (imdsinfo->coords) { + msglen += 3*4*num_coords+IMDHEADERSIZE; + } + if (imdsinfo->velocities) { + msglen += 3*4*num_coords+IMDHEADERSIZE; + } + if (imdsinfo->forces) { + msglen += 3*4*num_coords+IMDHEADERSIZE; + } + msgdata = new char[msglen]; + if (me == 0) { /* set up incoming socket on MPI rank 0. */ imdsock_init(); @@ -512,7 +604,7 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : if (imd_terminate) error->all(FLERR,"LAMMPS Terminated on error in IMD."); - /* storage required to communicate a single coordinate or force. */ + /* storage required to communicate a single coordinate, velocity, or force. */ size_one = sizeof(struct commdata); #if defined(LAMMPS_ASYNC_IMD) @@ -546,7 +638,6 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : *********************************/ FixIMD::~FixIMD() { - #if defined(LAMMPS_ASYNC_IMD) if (me == 0) { pthread_mutex_lock(&write_mutex); @@ -561,13 +652,17 @@ FixIMD::~FixIMD() pthread_cond_destroy(&write_cond); } #endif - auto hashtable = (taginthash_t *)idmap; - memory->destroy(comm_buf); - memory->destroy(force_buf); + memory->destroy(coord_data); + memory->destroy(vel_data); + memory->destroy(force_data); + + memory->destroy(msgdata); + memory->destroy(recv_force_buf); taginthash_destroy(hashtable); delete hashtable; free(rev_idmap); + free(imdsinfo); // close sockets imdsock_shutdown(clientsock); imdsock_destroy(clientsock); @@ -583,6 +678,7 @@ int FixIMD::setmask() int mask = 0; mask |= POST_FORCE; mask |= POST_FORCE_RESPA; + mask |= END_OF_STEP; return mask; } @@ -611,6 +707,7 @@ int FixIMD::reconnect() } else { fprintf(screen,"Waiting for IMD connection on port %d. Transfer rate %d.\n",imd_port, imd_trate); } + fflush(screen); } connect_msg = 0; clientsock = nullptr; @@ -637,7 +734,8 @@ int FixIMD::reconnect() return 0; } else { /* check endianness and IMD protocol version. */ - if (imd_handshake(clientsock)) { + if ((imd_version == 2 && imd_handshake_v2(clientsock)) || + (imd_version == 3 && imd_handshake_v3(clientsock, imdsinfo))) { if (screen) fprintf(screen, "IMD handshake error. Dropping connection.\n"); imdsock_destroy(clientsock); @@ -680,9 +778,21 @@ void FixIMD::setup(int) if (mask[i] & groupbit) ++nme; MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); - memory->destroy(comm_buf); + maxbuf = nmax*size_one; - comm_buf = (void *) memory->smalloc(maxbuf,"imd:comm_buf"); + + if (imdsinfo->coords) { + memory->destroy(coord_data); + coord_data = (void *) memory->smalloc(maxbuf,"imd:coord_data"); + } + if (imdsinfo->velocities) { + memory->destroy(vel_data); + vel_data = (void *) memory->smalloc(maxbuf,"imd:vel_data"); + } + if (imdsinfo->forces) { + memory->destroy(force_data); + force_data = (void *) memory->smalloc(maxbuf,"imd:force_data"); + } connect_msg = 1; reconnect(); @@ -697,11 +807,11 @@ void FixIMD::setup(int) idmap = (void *)hashtable; int tmp, ndata; - auto buf = static_cast(comm_buf); + auto buf = static_cast(coord_data); if (me == 0) { - MPI_Status status; - MPI_Request request; + std::vector statuses; + std::vector requests; auto taglist = new tagint[num_coords]; int numtag=0; /* counter to map atom tags to a 0-based consecutive index list */ @@ -714,15 +824,39 @@ void FixIMD::setup(int) /* loop over procs to receive remote data */ for (i=1; i < comm->nprocs; ++i) { - MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); - MPI_Send(&tmp, 0, MPI_INT, i, 0, world); - MPI_Wait(&request, &status); - MPI_Get_count(&status, MPI_BYTE, &ndata); - ndata /= size_one; + /* We're assuming tags are consistent across x,v,f */ + bool tag_recvd = false; + statuses.clear(); + requests.clear(); - for (j=0; j < ndata; ++j) { - taglist[numtag] = buf[j].tag; - ++numtag; + if (imdsinfo->coords) { + requests.push_back(MPI_Request()); + MPI_Irecv(coord_data, maxbuf, MPI_BYTE, i, 0, world, &requests.back()); + } + if (imdsinfo->velocities) { + requests.push_back(MPI_Request()); + MPI_Irecv(vel_data, maxbuf, MPI_BYTE, i, 0, world, &requests.back()); + } + if (imdsinfo->forces) { + requests.push_back(MPI_Request()); + MPI_Irecv(vel_data, maxbuf, MPI_BYTE, i, 0, world, &requests.back()); + } + statuses.resize(requests.size()); + MPI_Send(&tmp, 0, MPI_INT, i, 0, world); + MPI_Waitall(requests.size(), requests.data(), statuses.data()); + + for (size_t k=0; k < statuses.size(); ++k) { + if (!tag_recvd) { + MPI_Get_count(&statuses[k], MPI_BYTE, &ndata); + ndata /= size_one; + for (j=0; j < ndata; ++j) { + taglist[numtag] = buf[j].tag; + ++numtag; + } + tag_recvd = true; + } else { + break; + } } } @@ -747,7 +881,15 @@ void FixIMD::setup(int) } /* blocking receive to wait until it is our turn to send data. */ MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, MPI_STATUS_IGNORE); - MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world); + if (imdsinfo->coords) { + MPI_Rsend(coord_data, nme*size_one, MPI_BYTE, 0, 0, world); + } + if (imdsinfo->velocities) { + MPI_Rsend(vel_data, nme*size_one, MPI_BYTE, 0, 0, world); + } + if (imdsinfo->forces) { + MPI_Rsend(force_data, nme*size_one, MPI_BYTE, 0, 0, world); + } } } @@ -795,6 +937,39 @@ void FixIMD::ioworker() * Send coodinates, energies, and add IMD forces to atoms. */ void FixIMD::post_force(int /*vflag*/) { + fflush(screen); + if (imd_version == 2) { + handle_step_v2(); + } + else if (imd_version == 3) { + handle_client_input_v3(); + } + + } + +/* ---------------------------------------------------------------------- */ +void FixIMD::post_force_respa(int vflag, int ilevel, int /*iloop*/) +{ + /* only process IMD on the outmost RESPA level. */ + if (ilevel == nlevels_respa-1) post_force(vflag); +} + +void FixIMD::end_of_step() +{ + if (imd_version == 3 && update->ntimestep % imd_trate == 0) { + handle_output_v3(); + } +} + +/* ---------------------------------------------------------------------- */ +/* local memory usage. approximately. */ +double FixIMD::memory_usage() +{ + return static_cast(num_coords+maxbuf+imd_forces)*size_one; +} + +void FixIMD::handle_step_v2() { + /* check for reconnect */ if (imd_inactive) { reconnect(); @@ -825,22 +1000,12 @@ void FixIMD::post_force(int /*vflag*/) switch(msg) { - case IMD_GO: - if (screen) - fprintf(screen, "Ignoring unexpected IMD_GO message.\n"); - break; - - case IMD_IOERROR: - if (screen) - fprintf(screen, "IMD connection lost.\n"); - /* fallthrough */ - case IMD_DISCONNECT: { /* disconnect from client. wait for new connection. */ imd_paused = 0; imd_forces = 0; - memory->destroy(force_buf); - force_buf = nullptr; + memory->destroy(recv_force_buf); + recv_force_buf = nullptr; imdsock_destroy(clientsock); clientsock = nullptr; if (screen) @@ -866,9 +1031,9 @@ void FixIMD::post_force(int /*vflag*/) case IMD_PAUSE: /* pause the running simulation. wait for second IMD_PAUSE to continue. */ if (imd_paused) { - if (screen) - fprintf(screen, "Continuing run on IMD client request.\n"); - imd_paused = 0; + if (screen) + fprintf(screen, "Continuing run on IMD client request.\n"); + imd_paused = 0; } else { if (screen) fprintf(screen, "Pausing run on IMD client request.\n"); @@ -884,31 +1049,18 @@ void FixIMD::post_force(int /*vflag*/) fprintf(screen, "IMD client requested change of transfer rate. Now it is %d.\n", imd_trate); break; - case IMD_ENERGIES: { - IMDEnergies dummy_energies; - imd_recv_energies(clientsock, &dummy_energies); - break; - } - - case IMD_FCOORDS: { - auto dummy_coords = new float[3*length]; - imd_recv_fcoords(clientsock, length, dummy_coords); - delete[] dummy_coords; - break; - } - case IMD_MDCOMM: { auto imd_tags = new int32[length]; auto imd_fdat = new float[3*length]; imd_recv_mdcomm(clientsock, length, imd_tags, imd_fdat); if (imd_forces < length) { /* grow holding space for forces, if needed. */ - memory->destroy(force_buf); - force_buf = (void *) memory->smalloc((bigint)length*size_one, - "imd:force_buf"); + memory->destroy(recv_force_buf); + recv_force_buf = (void *) memory->smalloc((bigint)length*size_one, + "imd:recv_force_buf"); } imd_forces = length; - buf = static_cast(force_buf); + buf = static_cast(recv_force_buf); /* compare data to hash table */ for (int ii=0; ii < length; ++ii) { @@ -943,12 +1095,12 @@ void FixIMD::post_force(int /*vflag*/) /* check if we need to readjust the forces comm buffer on the receiving nodes. */ if (me != 0) { if (old_imd_forces < imd_forces) { /* grow holding space for forces, if needed. */ - if (force_buf != nullptr) - memory->sfree(force_buf); - force_buf = memory->smalloc((bigint)imd_forces*size_one, "imd:force_buf"); + if (recv_force_buf != nullptr) + memory->sfree(recv_force_buf); + recv_force_buf = memory->smalloc((bigint)imd_forces*size_one, "imd:recv_force_buf"); } } - MPI_Bcast(force_buf, imd_forces*size_one, MPI_BYTE, 0, world); + MPI_Bcast(recv_force_buf, imd_forces*size_one, MPI_BYTE, 0, world); } /* Check if we need to communicate coordinates to the client. @@ -961,7 +1113,7 @@ void FixIMD::post_force(int /*vflag*/) if (update->ntimestep % imd_trate) { if (imd_forces > 0) { double **f = atom->f; - buf = static_cast(force_buf); + buf = static_cast(recv_force_buf); /* XXX. this is in principle O(N**2) == not good. * however we assume for now that the number of atoms @@ -989,27 +1141,25 @@ void FixIMD::post_force(int /*vflag*/) MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); if (nmax*size_one > maxbuf) { - memory->destroy(comm_buf); + memory->destroy(coord_data); maxbuf = nmax*size_one; - comm_buf = memory->smalloc(maxbuf,"imd:comm_buf"); + coord_data = memory->smalloc(maxbuf,"imd:coord_data"); } int tmp, ndata; - buf = static_cast(comm_buf); + buf = static_cast(coord_data); if (me == 0) { MPI_Status status; MPI_Request request; /* collect data into new array. we bypass the IMD API to save * us one extra copy of the data. */ - msglen = 3*sizeof(float)*num_coords+IMDHEADERSIZE; - msgdata = new char[msglen]; imd_fill_header((IMDheader *)msgdata, IMD_FCOORDS, num_coords); /* array pointer, to the offset where we receive the coordinates. */ - auto recvcoord = (float *) (msgdata+IMDHEADERSIZE); + auto recvcoord = (float *) (msgdata+IMDHEADERSIZE); /* add local data */ - if (unwrap_flag) { + if (imdsinfo->unwrap) { double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; @@ -1049,10 +1199,9 @@ void FixIMD::post_force(int /*vflag*/) } } } - /* loop over procs to receive remote data */ for (i=1; i < comm->nprocs; ++i) { - MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); + MPI_Irecv(coord_data, maxbuf, MPI_BYTE, i, 0, world, &request); MPI_Send(&tmp, 0, MPI_INT, i, 0, world); MPI_Wait(&request, &status); MPI_Get_count(&status, MPI_BYTE, &ndata); @@ -1067,7 +1216,6 @@ void FixIMD::post_force(int /*vflag*/) } } } - /* done collecting frame data now communicate with IMD client. */ #if defined(LAMMPS_ASYNC_IMD) @@ -1081,13 +1229,12 @@ void FixIMD::post_force(int /*vflag*/) if (clientsock && imdsock_selwrite(clientsock,0)) { imd_writen(clientsock, msgdata, msglen); } - delete[] msgdata; #endif } else { /* copy coordinate data into communication buffer */ nme = 0; - if (unwrap_flag) { + if (imdsinfo->unwrap) { double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; @@ -1128,23 +1275,439 @@ void FixIMD::post_force(int /*vflag*/) } /* blocking receive to wait until it is our turn to send data. */ MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, MPI_STATUS_IGNORE); - MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world); + MPI_Rsend(coord_data, nme*size_one, MPI_BYTE, 0, 0, world); } - } - -/* ---------------------------------------------------------------------- */ -void FixIMD::post_force_respa(int vflag, int ilevel, int /*iloop*/) -{ - /* only process IMD on the outmost RESPA level. */ - if (ilevel == nlevels_respa-1) post_force(vflag); } -/* ---------------------------------------------------------------------- */ -/* local memory usage. approximately. */ -double FixIMD::memory_usage() -{ - return static_cast(num_coords+maxbuf+imd_forces)*size_one; +void FixIMD::handle_client_input_v3() { + // IMDV3 + + /* check for reconnect */ + if (imd_inactive) { + reconnect(); + MPI_Bcast(&imd_inactive, 1, MPI_INT, 0, world); + MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); + if (imd_terminate) + error->all(FLERR,"LAMMPS terminated on error in setting up IMD connection."); + if (imd_inactive) + return; /* IMD client has detached and not yet come back. do nothing. */ + } + + struct commdata *buf; + int nlocal = atom->nlocal; + int *mask = atom->mask; + tagint *tag = atom->tag; + double **f = atom->f; + + if (me == 0) { + /* process all pending incoming data. */ + int imd_paused=0; + while ((imdsock_selread(clientsock, 0) > 0) || imd_paused) { + /* if something requested to turn off IMD while paused get out */ + if (imd_inactive) break; + + int32 length; + int msg = imd_recv_header(clientsock, &length); + + switch(msg) { + + case IMD_DISCONNECT: { + /* disconnect from client. wait for new connection. */ + imd_paused = 0; + imd_forces = 0; + memory->destroy(recv_force_buf); + recv_force_buf = nullptr; + imdsock_destroy(clientsock); + clientsock = nullptr; + if (screen) + fprintf(screen, "IMD client detached. LAMMPS run continues.\n"); + + connect_msg = 1; + reconnect(); + if (imd_terminate) imd_inactive = 1; + break; + } + + case IMD_KILL: + /* stop the simulation job and shutdown IMD */ + if (screen) + fprintf(screen, "IMD client requested termination of run.\n"); + imd_inactive = 1; + imd_terminate = 1; + imd_paused = 0; + imdsock_destroy(clientsock); + clientsock = nullptr; + break; + + case IMD_PAUSE: + /* pause the running simulation. wait for second IMD_PAUSE to continue. */ + if (!imd_paused) { + if (screen) + fprintf(screen, "Pausing run on IMD client request.\n"); + imd_paused = 1; + } else { + // Pause in IMDv3 is idempotent + continue; + } + break; + + case IMD_RESUME: + /* resume the running simulation. */ + if (imd_paused) { + if (screen) + fprintf(screen, "Continuing run on IMD client request.\n"); + imd_paused = 0; + } else { + // Resume in IMDv3 is idempotent + continue; + } + break; + + case IMD_TRATE: + /* change the IMD transmission data rate */ + if (length > 0) + imd_trate = length; + if (screen) + fprintf(screen, "IMD client requested change of transfer rate. Now it is %d.\n", imd_trate); + break; + + case IMD_MDCOMM: { + auto imd_tags = new int32[length]; + auto imd_fdat = new float[3*length]; + imd_recv_mdcomm(clientsock, length, imd_tags, imd_fdat); + + if (imd_forces < length) { /* grow holding space for forces, if needed. */ + memory->destroy(recv_force_buf); + recv_force_buf = (void *) memory->smalloc((bigint)length*size_one, + "imd:recv_force_buf"); + } + imd_forces = length; + buf = static_cast(recv_force_buf); + + /* compare data to hash table */ + for (int ii=0; ii < length; ++ii) { + buf[ii].tag = rev_idmap[imd_tags[ii]]; + buf[ii].x = imd_fdat[3*ii]; + buf[ii].y = imd_fdat[3*ii+1]; + buf[ii].z = imd_fdat[3*ii+2]; + } + delete[] imd_tags; + delete[] imd_fdat; + break; + } + case IMD_WAIT: { + /* Change IMD waiting behavior mid-session */ + if (length) { + nowait_flag = 0; + } + else { + nowait_flag = 1; + } + break; + } + + default: + if (screen) + fprintf(screen, "Unhandled incoming IMD message #%d. length=%d\n", msg, length); + break; + } + } + } + + /* update all tasks with current settings. */ + int old_imd_forces = imd_forces; + MPI_Bcast(&imd_trate, 1, MPI_INT, 0, world); + MPI_Bcast(&imd_inactive, 1, MPI_INT, 0, world); + MPI_Bcast(&imd_forces, 1, MPI_INT, 0, world); + MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); + if (imd_terminate) + error->all(FLERR,"LAMMPS terminated on IMD request."); + + if (imd_forces > 0) { + /* check if we need to readjust the forces comm buffer on the receiving nodes. */ + if (me != 0) { + if (old_imd_forces < imd_forces) { /* grow holding space for forces, if needed. */ + if (recv_force_buf != nullptr) + memory->sfree(recv_force_buf); + recv_force_buf = memory->smalloc((bigint)imd_forces*size_one, "imd:recv_force_buf"); + } + } + MPI_Bcast(recv_force_buf, imd_forces*size_one, MPI_BYTE, 0, world); + } + + /* Check if we need to communicate coordinates to the client. + * Tuning imd_trate allows to keep the overhead for IMD low + * at the expense of a more jumpy display. Rather than using + * end_of_step() we do everything here in one go. + * + * If we don't communicate, only check if we have forces + * stored away and apply them. */ + if (imd_forces > 0) { + buf = static_cast(recv_force_buf); + + /* XXX. this is in principle O(N**2) == not good. + * however we assume for now that the number of atoms + * that we manipulate via IMD will be small compared + * to the total system size, so we don't hurt too much. */ + for (int j=0; j < imd_forces; ++j) { + for (int i=0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + if (buf[j].tag == tag[i]) { + f[i][0] += imd_fscale*buf[j].x; + f[i][1] += imd_fscale*buf[j].y; + f[i][2] += imd_fscale*buf[j].z; + } + } + } + } + } +} + +void FixIMD::handle_output_v3() { + + tagint *tag = atom->tag; + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + imageint *image = atom->image; + int nlocal = atom->nlocal; + int *mask = atom->mask; + struct commdata *buf; + + // Only main process should use: + float *global_coords = nullptr; + float *global_vel = nullptr; + float *global_force = nullptr; + + // Prepare offsets in outgoing buffer + // Fill what we can wihtout collecting from other processes + if (me == 0) { + int offset = 0; + if (imdsinfo->time) { + imd_fill_header((IMDheader *)msgdata, IMD_TIME, 1); + double dt = update->dt; + + double currtime = update->atime + ((update->ntimestep - update->atimestep) * update->dt); + long long currstep = update->ntimestep; + char *time = (msgdata+IMDHEADERSIZE); + + memcpy(time, &dt, sizeof(dt)); + memcpy(time+sizeof(dt), &currtime, sizeof(currtime)); + memcpy(time+sizeof(dt)+sizeof(currtime), &currstep, sizeof(currstep)); + offset += IMDHEADERSIZE+sizeof(dt)+sizeof(currtime)+sizeof(currstep); + } + if (imdsinfo->box) { + imd_fill_header((IMDheader *)(msgdata + offset), IMD_BOX, 1); + // Get triclinic box vectors + float *box = (float *)(msgdata+offset+IMDHEADERSIZE); + box[0] = domain->h[0]; + box[1] = 0.0; + box[2] = 0.0; + box[3] = domain->h[5]; + box[4] = domain->h[1]; + box[5] = 0.0; + box[6] = domain->h[4]; + box[7] = domain->h[3]; + box[8] = domain->h[2]; + + offset += (9*4)+IMDHEADERSIZE; + + } + if (imdsinfo->coords) { + imd_fill_header((IMDheader *)(msgdata + offset), IMD_FCOORDS, num_coords); + global_coords = (float *) (msgdata + offset + IMDHEADERSIZE); + offset += 3*4*num_coords+IMDHEADERSIZE; + } + if (imdsinfo->velocities) { + imd_fill_header((IMDheader *)(msgdata + offset), IMD_VELOCITIES, num_coords); + global_vel = (float *) (msgdata + offset + IMDHEADERSIZE); + offset += 3*4*num_coords+IMDHEADERSIZE; + } + if (imdsinfo->forces) { + imd_fill_header((IMDheader *)(msgdata + offset), IMD_FORCES, num_coords); + global_force = (float *) (msgdata + offset + IMDHEADERSIZE); + offset += 3*4*num_coords+IMDHEADERSIZE; + } + } + + int ntotal, nmax, nme=0; + for (int i=0; i < nlocal; ++i) + if (mask[i] & groupbit) ++nme; + + // Atoms per proc + int *recvcounts = nullptr; + // Displacements in recv for each proc + int *displs = nullptr; + + + if (me == 0) { + recvcounts = new int[comm->nprocs]; + displs = new int[comm->nprocs]; + } + + MPI_Gather(&nme, 1, MPI_INT, recvcounts, 1, MPI_INT, 0, world); + + if (me == 0) { + displs[0] = 0; + ntotal = recvcounts[0]; + for (int i = 1; i < comm->nprocs; ++i) { + displs[i] = displs[i - 1] + recvcounts[i - 1]; + ntotal += recvcounts[i]; + } + } + + if (imdsinfo->coords) { + commdata *recvcoord = nullptr; + memory->destroy(coord_data); + coord_data = memory->smalloc(nme*size_one, "imd:coord_data"); + buf = static_cast(coord_data); + int idx = 0; + if (imdsinfo->unwrap) { + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + double xy = domain->xy; + double xz = domain->xz; + double yz = domain->yz; + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + int ix = (image[i] & IMGMASK) - IMGMAX; + int iy = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; + int iz = (image[i] >> IMG2BITS) - IMGMAX; + + if (domain->triclinic) { + buf[idx].tag = tag[i]; + buf[idx].x = x[i][0]; + ix * xprd + iy * xy + iz * xz; + buf[idx].y = x[i][1]; + iy * yprd + iz * yz; + buf[idx].z = x[i][2]; + iz * zprd; + } + else { + buf[idx].tag = tag[i]; + buf[idx].x = x[i][0]; + ix * xprd; + buf[idx].y = x[i][1]; + iy * yprd; + buf[idx].z = x[i][2]; + iz * zprd; + } + ++idx; + } + } + } + else { + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + buf[idx].tag = tag[i]; + buf[idx].x = x[i][0]; + buf[idx].y = x[i][1]; + buf[idx].z = x[i][2]; + ++idx; + } + } + } + if (me == 0) { + recvcoord = new commdata[ntotal]; + } + MPI_Gatherv(buf, nme, MPI_CommData, + recvcoord, recvcounts, displs, MPI_CommData, 0, world); + if (me == 0) { + // Sort the coordinates by tag, place in global_coords + for (int i = 0; i < comm->nprocs; ++i) { + for (int j = 0; j < recvcounts[i]; ++j) { + int idx = displs[i]+j; + const tagint t = 3*taginthash_lookup((taginthash_t *)idmap, recvcoord[idx].tag); + if (t != 3*HASH_FAIL) { + global_coords[t] = recvcoord[idx].x; + global_coords[t + 1] = recvcoord[idx].y; + global_coords[t + 2] = recvcoord[idx].z; + } + } + } + } + } + if (imdsinfo->velocities) { + commdata *recvvel = nullptr; + memory->destroy(vel_data); + vel_data = memory->smalloc(nme*size_one, "imd:vel_data"); + buf = static_cast(vel_data); + int idx = 0; + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + buf[idx].tag = tag[i]; + buf[idx].x = v[i][0]; + buf[idx].y = v[i][1]; + buf[idx].z = v[i][2]; + ++idx; + } + } + if (me == 0) { + recvvel = new commdata[ntotal]; + } + MPI_Gatherv(buf, nme, MPI_CommData, + recvvel, recvcounts, displs, MPI_CommData, 0, world); + if (me == 0) { + // Sort the coordinates by tag, place in global_vels + for (int i = 0; i < comm->nprocs; ++i) { + for (int j = 0; j < recvcounts[i]; ++j) { + int idx = displs[i]+j; + const tagint t = 3*taginthash_lookup((taginthash_t *)idmap, recvvel[idx].tag); + if (t != 3*HASH_FAIL) { + global_vel[t] = recvvel[idx].x; + global_vel[t + 1] = recvvel[idx].y; + global_vel[t + 2] = recvvel[idx].z; + } + } + } + } + } + if (imdsinfo->forces) { + commdata *recvforce = nullptr; + memory->destroy(force_data); + force_data = memory->smalloc(nme*size_one, "imd:force_data"); + buf = static_cast(force_data); + int idx = 0; + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + buf[idx].tag = tag[i]; + buf[idx].x = f[i][0]; + buf[idx].y = f[i][1]; + buf[idx].z = f[i][2]; + ++idx; + } + } + if (me == 0) { + recvforce = new commdata[ntotal]; + } + MPI_Gatherv(buf, nme, MPI_CommData, + recvforce, recvcounts, displs, MPI_CommData, 0, world); + if (me == 0) { + // Sort the coordinates by tag, place in global_coords + for (int i = 0; i < comm->nprocs; ++i) { + for (int j = 0; j < recvcounts[i]; ++j) { + int idx = displs[i]+j; + const tagint t = 3*taginthash_lookup((taginthash_t *)idmap, recvforce[idx].tag); + if (t != 3*HASH_FAIL) { + global_force[t] = recvforce[idx].x; + global_force[t + 1] = recvforce[idx].y; + global_force[t + 2] = recvforce[idx].z; + } + } + } + } + } + +/* done collecting frame data now communicate with IMD client. */ + +#if defined(LAMMPS_ASYNC_IMD) + /* wake up i/o worker thread and release lock on i/o buffer + * we can go back to our MD and let the i/o thread do the rest */ + buf_has_data=1; + pthread_cond_signal(&write_cond); + pthread_mutex_unlock(&write_mutex); +#else + /* send coordinate data, if client is able to accept */ + if (clientsock && imdsock_selwrite(clientsock,0)) { + imd_writen(clientsock, msgdata, msglen); + } +#endif } /* End of FixIMD class implementation. */ @@ -1401,13 +1964,35 @@ static int32 imd_writen(void *s, const char *ptr, int32 n) { return n; } -int imd_handshake(void *s) { +int imd_handshake_v2(void *s) { IMDheader header; imd_fill_header(&header, IMD_HANDSHAKE, 1); - header.length = IMDVERSION; /* Not byteswapped! */ + header.length = 2; /* Not byteswapped! */ return (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE); } +int imd_handshake_v3(void *s, IMDSessionInfo *imdsinfo) { + IMDheader header; + imd_fill_header(&header, IMD_HANDSHAKE, 1); + header.length = 3; /* Not byteswapped so client can determine native endinaness */ + + if (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE) return -1; + + imd_fill_header(&header, IMD_SESSIONINFO, 7); + unsigned char body[7] = {0}; + body[0] = imdsinfo->time; + body[1] = imdsinfo->energies; + body[2] = imdsinfo->box; + body[3] = imdsinfo->coords; + body[4] = !imdsinfo->unwrap; + body[5] = imdsinfo->velocities; + body[6] = imdsinfo->forces; + + if (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE || + imd_writen(s, (char *)&body, 7) != 7) return -1; + return 0; +} + /* The IMD receive functions */ IMDType imd_recv_header(void *s, int32 *length) { diff --git a/src/MISC/fix_imd.h b/src/MISC/fix_imd.h index 2312540e66..6b8778dbf0 100644 --- a/src/MISC/fix_imd.h +++ b/src/MISC/fix_imd.h @@ -56,6 +56,9 @@ FixStyle(imd,FixIMD); #include #endif +/* IMDv3 session information */ +struct IMDSessionInfo; + /* prototype for c wrapper that calls the real worker */ extern "C" void *fix_imd_ioworker(void *); @@ -69,8 +72,11 @@ class FixIMD : public Fix { void init() override; void setup(int) override; void post_force(int) override; + void end_of_step() override; void post_force_respa(int, int, int) override; double memory_usage() override; + // Fix nevery at 1, use trate to skip in 'end_of_step` + int nevery = 1; protected: int imd_port; @@ -80,13 +86,17 @@ class FixIMD : public Fix { int num_coords; // total number of atoms controlled by this fix int size_one; // bytes per atom in communication buffer. int maxbuf; // size of atom communication buffer. - void *comm_buf; // communication buffer + void *coord_data; // communication buffer for coordinates + void *vel_data; // communication buffer for velocities + void *force_data; // communication buffer for forces void *idmap; // hash for mapping atom indices to consistent order. tagint *rev_idmap; // list of the hash keys for reverse mapping. - int imd_forces; // number of forces communicated via IMD. - void *force_buf; // force data buffer - double imd_fscale; // scale factor for forces. in case VMD's units are off. + int imd_version; // version of the IMD protocol to be used. + + int imd_forces; // number of forces communicated via IMD. + void *recv_force_buf; // force data buffer + double imd_fscale; // scale factor for forces. in case VMD's units are off. int imd_inactive; // true if IMD connection stopped. int imd_terminate; // true if IMD requests termination of run. @@ -96,12 +106,20 @@ class FixIMD : public Fix { int nowait_flag; // true if LAMMPS should not wait with the execution for VMD. int connect_msg; // flag to indicate whether a "listen for connection message" is needed. + /* IMDv3-only */ + IMDSessionInfo *imdsinfo; // session information for IMDv3 + int me; // my MPI rank in this "world". int nlevels_respa; // flag to determine respa levels. int msglen; char *msgdata; + private: + void handle_step_v2(); + void handle_client_input_v3(); + void handle_output_v3(); + #if defined(LAMMPS_ASYNC_IMD) int buf_has_data; // flag to indicate to the i/o thread what to do. pthread_mutex_t write_mutex; // mutex for sending coordinates to i/o thread From 940308ba5918da2d6c0f29a9786630387461d715 Mon Sep 17 00:00:00 2001 From: ljwoods2 Date: Tue, 12 Nov 2024 14:38:41 -0700 Subject: [PATCH 02/37] run CI --- .github/workflows/check-cpp23.yml | 1 - .github/workflows/check-vla.yml | 1 - .github/workflows/codeql-analysis.yml | 1 - .github/workflows/compile-msvc.yml | 1 - .github/workflows/coverity.yml | 1 - .github/workflows/full-regression.yml | 1 - .github/workflows/kokkos-regression.yaml | 1 - .github/workflows/quick-regression.yml | 1 - .github/workflows/style-check.yml | 1 - .github/workflows/unittest-linux.yml | 1 - .github/workflows/unittest-macos.yml | 1 - 11 files changed, 11 deletions(-) diff --git a/.github/workflows/check-cpp23.yml b/.github/workflows/check-cpp23.yml index 2cd53f2208..5022b44e0a 100644 --- a/.github/workflows/check-cpp23.yml +++ b/.github/workflows/check-cpp23.yml @@ -14,7 +14,6 @@ on: jobs: build: name: Build with C++23 support enabled - if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/check-vla.yml b/.github/workflows/check-vla.yml index ab89018a3d..af1e269333 100644 --- a/.github/workflows/check-vla.yml +++ b/.github/workflows/check-vla.yml @@ -14,7 +14,6 @@ on: jobs: build: name: Build with -Werror=vla - if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/codeql-analysis.yml b/.github/workflows/codeql-analysis.yml index c7dd945f5f..47fdae311f 100644 --- a/.github/workflows/codeql-analysis.yml +++ b/.github/workflows/codeql-analysis.yml @@ -10,7 +10,6 @@ on: jobs: analyze: name: Analyze - if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest permissions: diff --git a/.github/workflows/compile-msvc.yml b/.github/workflows/compile-msvc.yml index 7560bc0549..4f04979925 100644 --- a/.github/workflows/compile-msvc.yml +++ b/.github/workflows/compile-msvc.yml @@ -18,7 +18,6 @@ concurrency: jobs: build: name: Windows Compilation Test - if: ${{ github.repository == 'lammps/lammps' }} runs-on: windows-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/coverity.yml b/.github/workflows/coverity.yml index c0c3e3f89a..7e38643b4c 100644 --- a/.github/workflows/coverity.yml +++ b/.github/workflows/coverity.yml @@ -9,7 +9,6 @@ on: jobs: analyze: name: Analyze - if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest container: image: lammps/buildenv:ubuntu20.04 diff --git a/.github/workflows/full-regression.yml b/.github/workflows/full-regression.yml index a6b5353b9b..683f65409d 100644 --- a/.github/workflows/full-regression.yml +++ b/.github/workflows/full-regression.yml @@ -12,7 +12,6 @@ jobs: build: name: Build LAMMPS # restrict to official LAMMPS repository - if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/kokkos-regression.yaml b/.github/workflows/kokkos-regression.yaml index 0756b080b0..6238f15c93 100644 --- a/.github/workflows/kokkos-regression.yaml +++ b/.github/workflows/kokkos-regression.yaml @@ -12,7 +12,6 @@ jobs: build: name: Build LAMMPS with Kokkos OpenMP # restrict to official LAMMPS repository - if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/quick-regression.yml b/.github/workflows/quick-regression.yml index 88794bfa0a..5325e4b6cb 100644 --- a/.github/workflows/quick-regression.yml +++ b/.github/workflows/quick-regression.yml @@ -16,7 +16,6 @@ jobs: build: name: Build LAMMPS # restrict to official LAMMPS repository - if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/style-check.yml b/.github/workflows/style-check.yml index 7be2c4fc46..e97163269c 100644 --- a/.github/workflows/style-check.yml +++ b/.github/workflows/style-check.yml @@ -18,7 +18,6 @@ concurrency: jobs: build: name: Programming Style Conformance - if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest steps: diff --git a/.github/workflows/unittest-linux.yml b/.github/workflows/unittest-linux.yml index ce98fcea35..04b1d1f328 100644 --- a/.github/workflows/unittest-linux.yml +++ b/.github/workflows/unittest-linux.yml @@ -18,7 +18,6 @@ concurrency: jobs: build: name: Linux Unit Test - if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/unittest-macos.yml b/.github/workflows/unittest-macos.yml index 0d478a9d6b..d3e242774c 100644 --- a/.github/workflows/unittest-macos.yml +++ b/.github/workflows/unittest-macos.yml @@ -18,7 +18,6 @@ concurrency: jobs: build: name: MacOS Unit Test - if: ${{ github.repository == 'lammps/lammps' }} runs-on: macos-13 env: CCACHE_DIR: ${{ github.workspace }}/.ccache From f7915109f97038e88a1ed35f2a801411339f55d7 Mon Sep 17 00:00:00 2001 From: ljwoods2 Date: Mon, 25 Nov 2024 15:22:02 -0700 Subject: [PATCH 03/37] minor bug fixes --- src/MISC/fix_imd.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/MISC/fix_imd.cpp b/src/MISC/fix_imd.cpp index c295fdfaf8..9192f8a24c 100644 --- a/src/MISC/fix_imd.cpp +++ b/src/MISC/fix_imd.cpp @@ -456,7 +456,7 @@ MPI_Datatype MPI_CommData; /*************************************************************** * create class and parse arguments in LAMMPS script. Syntax: - * fix ID group-ID imd [version (2|3)] [unwrap (on|off)] [fscale ] [time (on|off)] [box (on|off)] [coordinates (on|off)] [velocities (on|off)] [forces (on|off)] + * fix ID group-ID imd [trate ] [version (2|3)] [unwrap (on|off)] [fscale ] [time (on|off)] [box (on|off)] [coordinates (on|off)] [velocities (on|off)] [forces (on|off)] ***************************************************************/ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) @@ -917,8 +917,10 @@ void FixIMD::ioworker() pthread_exit(nullptr); } else if (buf_has_data > 0) { /* send coordinate data, if client is able to accept */ - if (clientsock && imdsock_selwrite(clientsock,0)) { - imd_writen(clientsock, msgdata, msglen); + if (imd_writen(clientsock, msgdata, msglen) != msglen) { + fprintf(screen,"Asynchronous I/O thread terminated on error in sending IMDFrame.\n"); + pthread_mutex_unlock(&write_mutex); + pthread_exit(nullptr); } delete[] msgdata; buf_has_data=0; @@ -1226,8 +1228,8 @@ void FixIMD::handle_step_v2() { pthread_mutex_unlock(&write_mutex); #else /* send coordinate data, if client is able to accept */ - if (clientsock && imdsock_selwrite(clientsock,0)) { - imd_writen(clientsock, msgdata, msglen); + if (imd_writen(clientsock, msgdata, msglen) != msglen) { + error->all(FLERR, "LAMMPS terminated on error in sending IMDFrame"); } #endif @@ -1703,9 +1705,9 @@ void FixIMD::handle_output_v3() { pthread_cond_signal(&write_cond); pthread_mutex_unlock(&write_mutex); #else - /* send coordinate data, if client is able to accept */ - if (clientsock && imdsock_selwrite(clientsock,0)) { - imd_writen(clientsock, msgdata, msglen); + /* send IMDFrame data, blocking until client accepts */ + if (imd_writen(clientsock, msgdata, msglen) != msglen) { + error->all(FLERR, "LAMMPS terminated on error in sending IMDFrame"); } #endif } From aeb9003890c866994a916ea89e278b7982ce7fbe Mon Sep 17 00:00:00 2001 From: ljwoods2 Date: Mon, 25 Nov 2024 18:22:44 -0700 Subject: [PATCH 04/37] bug fixes --- src/MISC/fix_imd.cpp | 218 +++++++++++++++++++++++++++++-------------- src/MISC/fix_imd.h | 2 + 2 files changed, 152 insertions(+), 68 deletions(-) diff --git a/src/MISC/fix_imd.cpp b/src/MISC/fix_imd.cpp index 9192f8a24c..6f2bb53a8c 100644 --- a/src/MISC/fix_imd.cpp +++ b/src/MISC/fix_imd.cpp @@ -76,7 +76,6 @@ negotiate an appropriate license for such distribution." #endif #include -#include using namespace LAMMPS_NS; using namespace FixConst; @@ -461,7 +460,6 @@ MPI_Datatype MPI_CommData; FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - if (narg < 4) error->all(FLERR,"Illegal fix imd command"); @@ -508,7 +506,9 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : vel_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); } else if (0 == strcmp(arg[iarg], "forces")) { force_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); - } else error->all(FLERR,"Unknown fix imd parameter"); + } else { + error->all(FLERR,"Unknown fix imd parameter"); + } ++iarg; ++iarg; } @@ -542,7 +542,6 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : imdsinfo->energies = false; } - bigint n = group->count(igroup); if (n > MAXSMALLINT) error->all(FLERR,"Too many atoms for fix imd"); num_coords = static_cast (n); @@ -568,25 +567,30 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : idmap = nullptr; rev_idmap = nullptr; + if (imd_version == 3) { + msglen = 0; + if (imdsinfo->time) { + msglen += 24+IMDHEADERSIZE; + } + if (imdsinfo->box) { + msglen += 9*4+IMDHEADERSIZE; + } + if (imdsinfo->coords) { + msglen += 3*4*num_coords+IMDHEADERSIZE; + } + if (imdsinfo->velocities) { + msglen += 3*4*num_coords+IMDHEADERSIZE; + } + if (imdsinfo->forces) { + msglen += 3*4*num_coords+IMDHEADERSIZE; + } + msgdata = new char[msglen]; + } + else { + msglen = 3*sizeof(float)*num_coords+IMDHEADERSIZE; + msgdata = new char[msglen]; + } - msglen = 0; - if (imdsinfo->time) { - msglen += 24+IMDHEADERSIZE; - } - if (imdsinfo->box) { - msglen += 9*4+IMDHEADERSIZE; - } - if (imdsinfo->coords) { - msglen += 3*4*num_coords+IMDHEADERSIZE; - } - if (imdsinfo->velocities) { - msglen += 3*4*num_coords+IMDHEADERSIZE; - } - if (imdsinfo->forces) { - msglen += 3*4*num_coords+IMDHEADERSIZE; - } - msgdata = new char[msglen]; - if (me == 0) { /* set up incoming socket on MPI rank 0. */ imdsock_init(); @@ -638,6 +642,7 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : *********************************/ FixIMD::~FixIMD() { + #if defined(LAMMPS_ASYNC_IMD) if (me == 0) { pthread_mutex_lock(&write_mutex); @@ -657,12 +662,12 @@ FixIMD::~FixIMD() memory->destroy(vel_data); memory->destroy(force_data); - memory->destroy(msgdata); + delete[] msgdata; memory->destroy(recv_force_buf); taginthash_destroy(hashtable); delete hashtable; free(rev_idmap); - free(imdsinfo); + delete imdsinfo; // close sockets imdsock_shutdown(clientsock); imdsock_destroy(clientsock); @@ -763,6 +768,103 @@ int FixIMD::reconnect() /* wait for IMD client (e.g. VMD) to respond, initialize communication * buffers and collect tag/id maps. */ void FixIMD::setup(int) +{ + if (imd_version == 2) { + setup_v2(); + } + else { + setup_v3(); + } +} + +void FixIMD::setup_v2() { + /* nme: number of atoms in group on this MPI task + * nmax: max number of atoms in group across all MPI tasks + * nlocal: all local atoms + */ + int i,j; + int nmax,nme,nlocal; + int *mask = atom->mask; + tagint *tag = atom->tag; + nlocal = atom->nlocal; + nme=0; + for (i=0; i < nlocal; ++i) + if (mask[i] & groupbit) ++nme; + + MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); + memory->destroy(coord_data); + maxbuf = nmax*size_one; + coord_data = (void *) memory->smalloc(maxbuf,"imd:coord_data"); + + connect_msg = 1; + reconnect(); + MPI_Bcast(&imd_inactive, 1, MPI_INT, 0, world); + MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); + if (imd_terminate) + error->all(FLERR,"LAMMPS terminated on error in setting up IMD connection."); + + /* initialize and build hashtable. */ + auto hashtable=new taginthash_t; + taginthash_init(hashtable, num_coords); + idmap = (void *)hashtable; + + int tmp, ndata; + auto buf = static_cast(coord_data); + + if (me == 0) { + MPI_Status status; + MPI_Request request; + auto taglist = new tagint[num_coords]; + int numtag=0; /* counter to map atom tags to a 0-based consecutive index list */ + + for (i=0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + taglist[numtag] = tag[i]; + ++numtag; + } + } + + /* loop over procs to receive remote data */ + for (i=1; i < comm->nprocs; ++i) { + MPI_Irecv(coord_data, maxbuf, MPI_BYTE, i, 0, world, &request); + MPI_Send(&tmp, 0, MPI_INT, i, 0, world); + MPI_Wait(&request, &status); + MPI_Get_count(&status, MPI_BYTE, &ndata); + ndata /= size_one; + + for (j=0; j < ndata; ++j) { + taglist[numtag] = buf[j].tag; + ++numtag; + } + } + + /* sort list of tags by value to have consistently the + * same list when running in parallel and build hash table. */ + id_sort(taglist, 0, num_coords-1); + for (i=0; i < num_coords; ++i) { + taginthash_insert(hashtable, taglist[i], i); + } + delete[] taglist; + + /* generate reverse index-to-tag map for communicating + * IMD forces back to the proper atoms */ + rev_idmap=taginthash_keys(hashtable); + } else { + nme=0; + for (i=0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + buf[nme].tag = tag[i]; + ++nme; + } + } + /* blocking receive to wait until it is our turn to send data. */ + MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, MPI_STATUS_IGNORE); + MPI_Rsend(coord_data, nme*size_one, MPI_BYTE, 0, 0, world); + } + + } + +void FixIMD::setup_v3() { /* nme: number of atoms in group on this MPI task * nmax: max number of atoms in group across all MPI tasks @@ -807,11 +909,24 @@ void FixIMD::setup(int) idmap = (void *)hashtable; int tmp, ndata; - auto buf = static_cast(coord_data); + + struct commdata *buf = nullptr; + if (imdsinfo->coords) { + buf = static_cast(coord_data); + } + else if (imdsinfo->velocities) { + buf = static_cast(vel_data); + } + else if (imdsinfo->forces) { + buf = static_cast(force_data); + } if (me == 0) { - std::vector statuses; - std::vector requests; + if (buf == nullptr) { + return; + } + MPI_Status status; + MPI_Request request; auto taglist = new tagint[num_coords]; int numtag=0; /* counter to map atom tags to a 0-based consecutive index list */ @@ -824,39 +939,15 @@ void FixIMD::setup(int) /* loop over procs to receive remote data */ for (i=1; i < comm->nprocs; ++i) { - /* We're assuming tags are consistent across x,v,f */ - bool tag_recvd = false; - statuses.clear(); - requests.clear(); - - if (imdsinfo->coords) { - requests.push_back(MPI_Request()); - MPI_Irecv(coord_data, maxbuf, MPI_BYTE, i, 0, world, &requests.back()); - } - if (imdsinfo->velocities) { - requests.push_back(MPI_Request()); - MPI_Irecv(vel_data, maxbuf, MPI_BYTE, i, 0, world, &requests.back()); - } - if (imdsinfo->forces) { - requests.push_back(MPI_Request()); - MPI_Irecv(vel_data, maxbuf, MPI_BYTE, i, 0, world, &requests.back()); - } - statuses.resize(requests.size()); + MPI_Irecv(coord_data, maxbuf, MPI_BYTE, i, 0, world, &request); MPI_Send(&tmp, 0, MPI_INT, i, 0, world); - MPI_Waitall(requests.size(), requests.data(), statuses.data()); + MPI_Wait(&request, &status); + MPI_Get_count(&status, MPI_BYTE, &ndata); + ndata /= size_one; - for (size_t k=0; k < statuses.size(); ++k) { - if (!tag_recvd) { - MPI_Get_count(&statuses[k], MPI_BYTE, &ndata); - ndata /= size_one; - for (j=0; j < ndata; ++j) { - taglist[numtag] = buf[j].tag; - ++numtag; - } - tag_recvd = true; - } else { - break; - } + for (j=0; j < ndata; ++j) { + taglist[numtag] = buf[j].tag; + ++numtag; } } @@ -881,19 +972,10 @@ void FixIMD::setup(int) } /* blocking receive to wait until it is our turn to send data. */ MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, MPI_STATUS_IGNORE); - if (imdsinfo->coords) { - MPI_Rsend(coord_data, nme*size_one, MPI_BYTE, 0, 0, world); - } - if (imdsinfo->velocities) { - MPI_Rsend(vel_data, nme*size_one, MPI_BYTE, 0, 0, world); - } - if (imdsinfo->forces) { - MPI_Rsend(force_data, nme*size_one, MPI_BYTE, 0, 0, world); - } + MPI_Rsend(coord_data, nme*size_one, MPI_BYTE, 0, 0, world); } } - /* worker threads for asynchronous i/o */ #if defined(LAMMPS_ASYNC_IMD) /* c bindings wrapper */ diff --git a/src/MISC/fix_imd.h b/src/MISC/fix_imd.h index 6b8778dbf0..03d242f32b 100644 --- a/src/MISC/fix_imd.h +++ b/src/MISC/fix_imd.h @@ -116,6 +116,8 @@ class FixIMD : public Fix { char *msgdata; private: + void setup_v2(); + void setup_v3(); void handle_step_v2(); void handle_client_input_v3(); void handle_output_v3(); From 3d2999194700bfb8d045c76785162ceabf38848f Mon Sep 17 00:00:00 2001 From: ljwoods2 Date: Tue, 17 Dec 2024 14:22:45 -0700 Subject: [PATCH 05/37] fix async IO bug --- src/MISC/fix_imd.cpp | 2 +- src/MISC/fix_imd.h | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/MISC/fix_imd.cpp b/src/MISC/fix_imd.cpp index 6f2bb53a8c..ea6200ab9e 100644 --- a/src/MISC/fix_imd.cpp +++ b/src/MISC/fix_imd.cpp @@ -1004,7 +1004,7 @@ void FixIMD::ioworker() pthread_mutex_unlock(&write_mutex); pthread_exit(nullptr); } - delete[] msgdata; + /* Don't delete msgdata- allow memory to persist. */ buf_has_data=0; pthread_mutex_unlock(&write_mutex); } else { diff --git a/src/MISC/fix_imd.h b/src/MISC/fix_imd.h index 03d242f32b..3d35376868 100644 --- a/src/MISC/fix_imd.h +++ b/src/MISC/fix_imd.h @@ -77,6 +77,7 @@ class FixIMD : public Fix { double memory_usage() override; // Fix nevery at 1, use trate to skip in 'end_of_step` int nevery = 1; + int imd_version; // version of the IMD protocol to be used. protected: int imd_port; @@ -92,8 +93,6 @@ class FixIMD : public Fix { void *idmap; // hash for mapping atom indices to consistent order. tagint *rev_idmap; // list of the hash keys for reverse mapping. - int imd_version; // version of the IMD protocol to be used. - int imd_forces; // number of forces communicated via IMD. void *recv_force_buf; // force data buffer double imd_fscale; // scale factor for forces. in case VMD's units are off. From 2c9ad698ca4f46084008c9825882925d2ed7a8d5 Mon Sep 17 00:00:00 2001 From: ljwoods2 Date: Tue, 17 Dec 2024 14:36:41 -0700 Subject: [PATCH 06/37] doc updates --- doc/src/fix_imd.rst | 36 ++++++++++++++++++++++++++++++++++-- 1 file changed, 34 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_imd.rst b/doc/src/fix_imd.rst index 520af505a1..d66da2c0d5 100644 --- a/doc/src/fix_imd.rst +++ b/doc/src/fix_imd.rst @@ -26,6 +26,24 @@ Syntax *nowait* arg = *on* or *off* off = LAMMPS waits to be connected to an IMD client before continuing (default) on = LAMMPS listens for an IMD client, but continues with the run + *version* arg = *2* or *3* + 2 = use IMD version 2 (default) + 3 = use IMD version 3. The subsequent keywords are only supported for version 3. + *time* arg = *on* or *off* + off = simulation time is not transmitted (default) + on = simulation time is transmitted. + *box* arg = *on* or *off* + off = simulation box data is not transmitted (default) + on = simulation box data is transmitted. + *coordinates* arg = *on* or *off* + off = atomic coordinates are not transmitted (default) + on = atomic coordinates are transmitted. + *velocities* arg = *on* or *off* + off = atomic velocities are not transmitted (default) + on = atomic velocities are transmitted. + *forces* arg = *on* or *off* + off = atomic forces are not transmitted (default) + on = atomic forces are transmitted. Examples """""""" @@ -41,10 +59,13 @@ Description This fix implements the "Interactive MD" (IMD) protocol which allows realtime visualization and manipulation of MD simulations through the IMD protocol, as initially implemented in VMD and NAMD. Specifically -it allows LAMMPS to connect an IMD client, for example the `VMD visualization program `_, so that it can monitor the progress of the +it allows LAMMPS to connect an IMD client, for example the `VMD visualization program `_ +(currently only supports IMDv2) or the +`Python IMDClient `_ (supports both IMDv2 and IMDv3), +so that it can monitor the progress of the simulation and interactively apply forces to selected atoms. -If LAMMPS is compiled with the pre-processor flag -DLAMMPS_ASYNC_IMD +If LAMMPS is compiled with the pre-processor flag -DLAMMPS_ASYNC_IMD (-DCMAKE_CXX_FLAGS="-DLAMMPS_ASYNC_IMD" in CMake), then fix imd will use POSIX threads to spawn a IMD communication thread on MPI rank 0 in order to offload data reading and writing from the main execution thread and potentially lower the inferred @@ -94,6 +115,13 @@ with different units or as a measure to tweak the forces generated by the manipulation of the IMD client, this option allows to make adjustments. +In `IMDv3 `_, the protocol has been extended to allow for the transmission +of simulation time, box dimensions, atomic coordinates, velocities, and +forces. The *version* keyword allows to select the version of the +protocol to be used. The *time*, *box*, *coordinates*, *velocities*, +and *forces* keywords allow to select which data is transmitted to the +IMD client. The default is to transmit all data. + To connect VMD to a listening LAMMPS simulation on the same machine with fix imd enabled, one needs to start VMD and load a coordinate or topology file that matches the fix group. When the VMD command @@ -129,6 +157,10 @@ screen output is active. .. _imdvmd: https://www.ks.uiuc.edu/Research/vmd/imd/ +.. _IMDClient: https://github.com/Becksteinlab/imdclient/tree/main/imdclient + +.. _IMDv3: https://imdclient.readthedocs.io/en/latest/protocol_v3.html + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" From 3c026df75d2a7a8f26c198d83cb6708b853db4bb Mon Sep 17 00:00:00 2001 From: ljwoods2 Date: Tue, 17 Dec 2024 14:36:59 -0700 Subject: [PATCH 07/37] formatting docs --- doc/src/fix_imd.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_imd.rst b/doc/src/fix_imd.rst index d66da2c0d5..a8400e49c8 100644 --- a/doc/src/fix_imd.rst +++ b/doc/src/fix_imd.rst @@ -41,7 +41,7 @@ Syntax *velocities* arg = *on* or *off* off = atomic velocities are not transmitted (default) on = atomic velocities are transmitted. - *forces* arg = *on* or *off* + *forces* arg = *on* or *off* off = atomic forces are not transmitted (default) on = atomic forces are transmitted. From cec242a421fb869303f5b75c6824c80e2efaf93a Mon Sep 17 00:00:00 2001 From: ljwoods2 Date: Mon, 6 Jan 2025 12:30:45 -0700 Subject: [PATCH 08/37] remove trailing whitespace --- doc/src/fix_imd.rst | 4 ++-- src/MISC/fix_imd.cpp | 37 ++++++++++++++++++------------------- 2 files changed, 20 insertions(+), 21 deletions(-) diff --git a/doc/src/fix_imd.rst b/doc/src/fix_imd.rst index a8400e49c8..96bcb52396 100644 --- a/doc/src/fix_imd.rst +++ b/doc/src/fix_imd.rst @@ -59,8 +59,8 @@ Description This fix implements the "Interactive MD" (IMD) protocol which allows realtime visualization and manipulation of MD simulations through the IMD protocol, as initially implemented in VMD and NAMD. Specifically -it allows LAMMPS to connect an IMD client, for example the `VMD visualization program `_ -(currently only supports IMDv2) or the +it allows LAMMPS to connect an IMD client, for example the `VMD visualization program `_ +(currently only supports IMDv2) or the `Python IMDClient `_ (supports both IMDv2 and IMDv3), so that it can monitor the progress of the simulation and interactively apply forces to selected atoms. diff --git a/src/MISC/fix_imd.cpp b/src/MISC/fix_imd.cpp index ea6200ab9e..5be6de2c2c 100644 --- a/src/MISC/fix_imd.cpp +++ b/src/MISC/fix_imd.cpp @@ -43,7 +43,7 @@ negotiate an appropriate license for such distribution." ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Axel Kohlmeyer (Temple U) + Contributing authors: Axel Kohlmeyer (Temple U), Lawson Woods (Arizona State U) IMD API, hash, and socket code written by: John E. Stone, Justin Gullingsrud, and James Phillips, (TCBG, Beckman Institute, UIUC) ------------------------------------------------------------------------- */ @@ -469,7 +469,7 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : /* default values for optional flags */ imd_version = 2; - + unwrap_flag = 0; nowait_flag = 0; connect_msg = 1; @@ -482,7 +482,7 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : int coord_flag = 1; int vel_flag = 1; int force_flag = 1; - + /* parse optional arguments */ int iarg = 4; while (iarg+1 < narg) { @@ -528,7 +528,7 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : imdsinfo->coords = true; imdsinfo->unwrap = unwrap_flag; imdsinfo->velocities = false; - imdsinfo->forces = false; + imdsinfo->forces = false; imdsinfo->energies = false; } @@ -551,7 +551,7 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : MPI_Type_commit(&MPI_CommData); MPI_Comm_rank(world,&me); - + /* initialize various imd state variables. */ clientsock = nullptr; localsock = nullptr; @@ -590,7 +590,7 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : msglen = 3*sizeof(float)*num_coords+IMDHEADERSIZE; msgdata = new char[msglen]; } - + if (me == 0) { /* set up incoming socket on MPI rank 0. */ imdsock_init(); @@ -864,7 +864,7 @@ void FixIMD::setup_v2() { } -void FixIMD::setup_v3() +void FixIMD::setup_v3() { /* nme: number of atoms in group on this MPI task * nmax: max number of atoms in group across all MPI tasks @@ -880,7 +880,7 @@ void FixIMD::setup_v3() if (mask[i] & groupbit) ++nme; MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); - + maxbuf = nmax*size_one; if (imdsinfo->coords) { @@ -1039,7 +1039,7 @@ void FixIMD::post_force_respa(int vflag, int ilevel, int /*iloop*/) } void FixIMD::end_of_step() -{ +{ if (imd_version == 3 && update->ntimestep % imd_trate == 0) { handle_output_v3(); } @@ -1240,7 +1240,7 @@ void FixIMD::handle_step_v2() { * us one extra copy of the data. */ imd_fill_header((IMDheader *)msgdata, IMD_FCOORDS, num_coords); /* array pointer, to the offset where we receive the coordinates. */ - auto recvcoord = (float *) (msgdata+IMDHEADERSIZE); + auto recvcoord = (float *) (msgdata+IMDHEADERSIZE); /* add local data */ if (imdsinfo->unwrap) { @@ -1435,7 +1435,7 @@ void FixIMD::handle_client_input_v3() { continue; } break; - + case IMD_RESUME: /* resume the running simulation. */ if (imd_paused) { @@ -1571,7 +1571,7 @@ void FixIMD::handle_output_v3() { if (imdsinfo->time) { imd_fill_header((IMDheader *)msgdata, IMD_TIME, 1); double dt = update->dt; - + double currtime = update->atime + ((update->ntimestep - update->atimestep) * update->dt); long long currstep = update->ntimestep; char *time = (msgdata+IMDHEADERSIZE); @@ -1594,9 +1594,9 @@ void FixIMD::handle_output_v3() { box[6] = domain->h[4]; box[7] = domain->h[3]; box[8] = domain->h[2]; - + offset += (9*4)+IMDHEADERSIZE; - + } if (imdsinfo->coords) { imd_fill_header((IMDheader *)(msgdata + offset), IMD_FCOORDS, num_coords); @@ -1631,12 +1631,11 @@ void FixIMD::handle_output_v3() { } MPI_Gather(&nme, 1, MPI_INT, recvcounts, 1, MPI_INT, 0, world); - if (me == 0) { displs[0] = 0; ntotal = recvcounts[0]; for (int i = 1; i < comm->nprocs; ++i) { - displs[i] = displs[i - 1] + recvcounts[i - 1]; + displs[i] = displs[i - 1] + recvcounts[i - 1]; ntotal += recvcounts[i]; } } @@ -1699,7 +1698,7 @@ void FixIMD::handle_output_v3() { int idx = displs[i]+j; const tagint t = 3*taginthash_lookup((taginthash_t *)idmap, recvcoord[idx].tag); if (t != 3*HASH_FAIL) { - global_coords[t] = recvcoord[idx].x; + global_coords[t] = recvcoord[idx].x; global_coords[t + 1] = recvcoord[idx].y; global_coords[t + 2] = recvcoord[idx].z; } @@ -2071,8 +2070,8 @@ int imd_handshake_v3(void *s, IMDSessionInfo *imdsinfo) { body[4] = !imdsinfo->unwrap; body[5] = imdsinfo->velocities; body[6] = imdsinfo->forces; - - if (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE || + + if (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE || imd_writen(s, (char *)&body, 7) != 7) return -1; return 0; } From 0dacedd8b0fddc2f70cc597e4289b188e9fcaee2 Mon Sep 17 00:00:00 2001 From: ljwoods2 Date: Wed, 8 Jan 2025 16:12:58 -0700 Subject: [PATCH 09/37] undo workflow changes --- .github/workflows/check-cpp23.yml | 1 + .github/workflows/check-vla.yml | 1 + .github/workflows/codeql-analysis.yml | 1 + .github/workflows/compile-msvc.yml | 1 + .github/workflows/coverity.yml | 1 + .github/workflows/full-regression.yml | 1 + .github/workflows/kokkos-regression.yaml | 1 + .github/workflows/quick-regression.yml | 1 + .github/workflows/style-check.yml | 1 + .github/workflows/unittest-linux.yml | 1 + .github/workflows/unittest-macos.yml | 1 + 11 files changed, 11 insertions(+) diff --git a/.github/workflows/check-cpp23.yml b/.github/workflows/check-cpp23.yml index 5022b44e0a..2cd53f2208 100644 --- a/.github/workflows/check-cpp23.yml +++ b/.github/workflows/check-cpp23.yml @@ -14,6 +14,7 @@ on: jobs: build: name: Build with C++23 support enabled + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/check-vla.yml b/.github/workflows/check-vla.yml index af1e269333..ab89018a3d 100644 --- a/.github/workflows/check-vla.yml +++ b/.github/workflows/check-vla.yml @@ -14,6 +14,7 @@ on: jobs: build: name: Build with -Werror=vla + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/codeql-analysis.yml b/.github/workflows/codeql-analysis.yml index 47fdae311f..c7dd945f5f 100644 --- a/.github/workflows/codeql-analysis.yml +++ b/.github/workflows/codeql-analysis.yml @@ -10,6 +10,7 @@ on: jobs: analyze: name: Analyze + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest permissions: diff --git a/.github/workflows/compile-msvc.yml b/.github/workflows/compile-msvc.yml index 4f04979925..7560bc0549 100644 --- a/.github/workflows/compile-msvc.yml +++ b/.github/workflows/compile-msvc.yml @@ -18,6 +18,7 @@ concurrency: jobs: build: name: Windows Compilation Test + if: ${{ github.repository == 'lammps/lammps' }} runs-on: windows-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/coverity.yml b/.github/workflows/coverity.yml index 7e38643b4c..c0c3e3f89a 100644 --- a/.github/workflows/coverity.yml +++ b/.github/workflows/coverity.yml @@ -9,6 +9,7 @@ on: jobs: analyze: name: Analyze + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest container: image: lammps/buildenv:ubuntu20.04 diff --git a/.github/workflows/full-regression.yml b/.github/workflows/full-regression.yml index 683f65409d..a6b5353b9b 100644 --- a/.github/workflows/full-regression.yml +++ b/.github/workflows/full-regression.yml @@ -12,6 +12,7 @@ jobs: build: name: Build LAMMPS # restrict to official LAMMPS repository + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/kokkos-regression.yaml b/.github/workflows/kokkos-regression.yaml index b3e685806b..a475c1eb7b 100644 --- a/.github/workflows/kokkos-regression.yaml +++ b/.github/workflows/kokkos-regression.yaml @@ -12,6 +12,7 @@ jobs: build: name: Build LAMMPS with Kokkos OpenMP # restrict to official LAMMPS repository + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/quick-regression.yml b/.github/workflows/quick-regression.yml index 5325e4b6cb..88794bfa0a 100644 --- a/.github/workflows/quick-regression.yml +++ b/.github/workflows/quick-regression.yml @@ -16,6 +16,7 @@ jobs: build: name: Build LAMMPS # restrict to official LAMMPS repository + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/style-check.yml b/.github/workflows/style-check.yml index e97163269c..7be2c4fc46 100644 --- a/.github/workflows/style-check.yml +++ b/.github/workflows/style-check.yml @@ -18,6 +18,7 @@ concurrency: jobs: build: name: Programming Style Conformance + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest steps: diff --git a/.github/workflows/unittest-linux.yml b/.github/workflows/unittest-linux.yml index 04b1d1f328..ce98fcea35 100644 --- a/.github/workflows/unittest-linux.yml +++ b/.github/workflows/unittest-linux.yml @@ -18,6 +18,7 @@ concurrency: jobs: build: name: Linux Unit Test + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest env: CCACHE_DIR: ${{ github.workspace }}/.ccache diff --git a/.github/workflows/unittest-macos.yml b/.github/workflows/unittest-macos.yml index d3e242774c..0d478a9d6b 100644 --- a/.github/workflows/unittest-macos.yml +++ b/.github/workflows/unittest-macos.yml @@ -18,6 +18,7 @@ concurrency: jobs: build: name: MacOS Unit Test + if: ${{ github.repository == 'lammps/lammps' }} runs-on: macos-13 env: CCACHE_DIR: ${{ github.workspace }}/.ccache From 02a8a9d706cf7675d507f99814340e9775d03bc6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 8 Jan 2025 20:30:48 -0500 Subject: [PATCH 10/37] make LAMMPS_ASYNC_IMD a CMake variable, too, and document it properly --- cmake/CMakeLists.txt | 10 ++++++++++ doc/src/Build_extras.rst | 34 +++++++++++++++++++++++++++++++++- doc/src/Build_package.rst | 1 + doc/src/fix_imd.rst | 24 ++++++++++++------------ 4 files changed, 56 insertions(+), 13 deletions(-) diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 38431e3bb5..e89a2a3a77 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -579,6 +579,16 @@ foreach(PKG_WITH_INCL KSPACE PYTHON ML-IAP VORONOI COLVARS ML-HDNNP MDI MOLFILE endif() endforeach() +# settings for misc packages and styles +if(PKG_MISC) + option(LAMMPS_ASYNC_IMD "Asynchronous IMD processing" OFF) + mark_as_advanced(LAMMPS_ASYNC_IMD) + if(LAMMPS_ASYNC_IMD) + target_compile_definitions(lammps PRIVATE -DLAMMPS_ASYNC_IMD) + message(STATUS "Using IMD in asynchronous mode") + endif() +endif() + # optionally enable building script wrappers using swig option(WITH_SWIG "Build scripting language wrappers with SWIG" OFF) if(WITH_SWIG) diff --git a/doc/src/Build_extras.rst b/doc/src/Build_extras.rst index 08e4cdbb7a..80b9da1bd5 100644 --- a/doc/src/Build_extras.rst +++ b/doc/src/Build_extras.rst @@ -48,6 +48,7 @@ This is the list of packages that may require additional steps. * :ref:`LEPTON ` * :ref:`MACHDYN ` * :ref:`MDI ` + * :ref:`MISC ` * :ref:`ML-HDNNP ` * :ref:`ML-IAP ` * :ref:`ML-PACE ` @@ -2019,7 +2020,7 @@ TBB and MKL. .. _mdi: MDI package ------------------------------ +----------- .. tabs:: @@ -2046,6 +2047,37 @@ MDI package ---------- +.. _misc: + +MISC package +------------ + +The :doc:`fix imd ` style in this package can be run either +synchronously (communication with IMD clients is done in the main +process) or asynchronously (the fix spawns a separate thread that can +communicate with IMD clients concurrently to the LAMMPS execution). + +.. tabs:: + + .. tab:: CMake build + + .. code-block:: bash + + -D LAMMPS_ASYNC_IMD=value # Run IMD server asynchronously + # value = no (default) or yes + + .. tab:: Traditional make + + To enable asynchronous mode the ``-DLAMMPS_ASYNC_IMD`` define + needs to be added to the ``LMP_INC`` variable in the + ``Makefile.machine`` you are using. For example: + + .. code-block:: make + + LMP_INC = -DLAMMPS_ASYNC_IMD -DLAMMPS_MEMALIGN=64 + +---------- + .. _molfile: MOLFILE package diff --git a/doc/src/Build_package.rst b/doc/src/Build_package.rst index 8b2da8ad30..c4c4889806 100644 --- a/doc/src/Build_package.rst +++ b/doc/src/Build_package.rst @@ -49,6 +49,7 @@ packages: * :ref:`LEPTON ` * :ref:`MACHDYN ` * :ref:`MDI ` + * :ref:`MISC ` * :ref:`ML-HDNNP ` * :ref:`ML-IAP ` * :ref:`ML-PACE ` diff --git a/doc/src/fix_imd.rst b/doc/src/fix_imd.rst index 1ec279a91b..158b44bdc6 100644 --- a/doc/src/fix_imd.rst +++ b/doc/src/fix_imd.rst @@ -58,19 +58,19 @@ Description This fix implements the "Interactive MD" (IMD) protocol which allows realtime visualization and manipulation of MD simulations through the -IMD protocol, as initially implemented in VMD and NAMD. Specifically -it allows LAMMPS to connect an IMD client, for example the `VMD visualization program `_ -(currently only supports IMDv2) or the -`Python IMDClient `_ (supports both IMDv2 and IMDv3), -so that it can monitor the progress of the -simulation and interactively apply forces to selected atoms. +IMD protocol, as initially implemented in VMD and NAMD. Specifically it +allows LAMMPS to connect an IMD client, for example the `VMD +visualization program `_ (currently only supports IMDv2) or the +`Python IMDClient `_ (supports both IMDv2 and IMDv3), so +that it can monitor the progress of the simulation and interactively +apply forces to selected atoms. -If LAMMPS is compiled with the pre-processor flag -DLAMMPS_ASYNC_IMD (-DCMAKE_CXX_FLAGS="-DLAMMPS_ASYNC_IMD" in CMake), -then fix imd will use POSIX threads to spawn a IMD communication -thread on MPI rank 0 in order to offload data reading and writing -from the main execution thread and potentially lower the inferred -latencies for slow communication links. This feature has only been -tested under linux. +If LAMMPS is compiled with the pre-processor flag +:ref:`-DLAMMPS_ASYNC_IMD <_misc>` then fix imd will use POSIX threads to +spawn an IMD communication thread on MPI rank 0 in order to offload data +exchange with the IMD client from the main execution thread and +potentially lower the inferred latencies for slow communication +links. This feature has only been tested under linux. The source code for this fix includes code developed by the Theoretical and Computational Biophysics Group in the Beckman Institute for Advanced From 6a363d441b9f6ed4fef39268719a459a01f15fd1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 8 Jan 2025 20:35:08 -0500 Subject: [PATCH 11/37] small formatting tweaks, add versionadded tag --- doc/src/fix_imd.rst | 35 +++++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/doc/src/fix_imd.rst b/doc/src/fix_imd.rst index 158b44bdc6..fd0fd2dc00 100644 --- a/doc/src/fix_imd.rst +++ b/doc/src/fix_imd.rst @@ -27,8 +27,13 @@ Syntax off = LAMMPS waits to be connected to an IMD client before continuing (default) on = LAMMPS listens for an IMD client, but continues with the run *version* arg = *2* or *3* - 2 = use IMD version 2 (default) - 3 = use IMD version 3. The subsequent keywords are only supported for version 3. + 2 = use IMD protocol version 2 (default) + 3 = use IMD protocol version 3. + + The following keywords are only supported for IMD protocol version 3. + + .. parsed-literal:: + *time* arg = *on* or *off* off = simulation time is not transmitted (default) on = simulation time is transmitted. @@ -115,12 +120,14 @@ with different units or as a measure to tweak the forces generated by the manipulation of the IMD client, this option allows to make adjustments. -In `IMDv3 `_, the protocol has been extended to allow for the transmission -of simulation time, box dimensions, atomic coordinates, velocities, and -forces. The *version* keyword allows to select the version of the -protocol to be used. The *time*, *box*, *coordinates*, *velocities*, -and *forces* keywords allow to select which data is transmitted to the -IMD client. The default is to transmit all data. +.. versionadded:: TBD + +In `IMDv3 `_, the IMD protocol has been extended to allow for +the transmission of simulation time, box dimensions, atomic coordinates, +velocities, and forces. The *version* keyword allows to select the +version of the protocol to be used. The *time*, *box*, *coordinates*, +*velocities*, and *forces* keywords allow to select which data is +transmitted to the IMD client. The default is to transmit all data. To connect VMD to a listening LAMMPS simulation on the same machine with fix imd enabled, one needs to start VMD and load a coordinate or @@ -179,14 +186,14 @@ This fix is part of the MISC package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` page for more info. -When used in combination with VMD, a topology or coordinate file has -to be loaded, which matches (in number and ordering of atoms) the -group the fix is applied to. The fix internally sorts atom IDs by -ascending integer value; in VMD (and thus the IMD protocol) those will -be assigned 0-based consecutive index numbers. +When used in combination with VMD, a topology or coordinate file has to +be loaded, which matches (in number and ordering of atoms) the group the +fix is applied to. The fix internally sorts atom IDs by ascending +integer value; in VMD (and thus the IMD protocol) those will be assigned +0-based consecutive index numbers. When using multiple active IMD connections at the same time, each -needs to use a different port number. +fix instance needs to use a different port number. Related commands """""""""""""""" From 1032c94c9f5c9e256dacf2ae73d7a33317d01179 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 8 Jan 2025 20:52:57 -0500 Subject: [PATCH 12/37] apply some more of LAMMPS' programming style guidelines and improve errors --- src/MISC/fix_imd.cpp | 38 ++++++++++++++++---------------------- 1 file changed, 16 insertions(+), 22 deletions(-) diff --git a/src/MISC/fix_imd.cpp b/src/MISC/fix_imd.cpp index 5be6de2c2c..493c3a4e0f 100644 --- a/src/MISC/fix_imd.cpp +++ b/src/MISC/fix_imd.cpp @@ -454,29 +454,31 @@ struct commdata { MPI_Datatype MPI_CommData; /*************************************************************** - * create class and parse arguments in LAMMPS script. Syntax: - * fix ID group-ID imd [trate ] [version (2|3)] [unwrap (on|off)] [fscale ] [time (on|off)] [box (on|off)] [coordinates (on|off)] [velocities (on|off)] [forces (on|off)] + * create class and parse arguments in LAMMPS script. ***************************************************************/ + FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) + Fix(lmp, narg, arg), localsock(nullptr), clientsock(nullptr), coord_data(nullptr), + vel_data(nullptr), force_data(nullptr), idmap(nullptr), rev_idmap(nullptr), + recv_force_buf(nullptr), imdsinfo(nullptr), msgdata(nullptr) { - if (narg < 4) - error->all(FLERR,"Illegal fix imd command"); + if (narg < 4) utils::missing_cmd_args(FLERR, "fix imd", error); imd_port = utils::inumeric(FLERR,arg[3],false,lmp); - if (imd_port < 1024) - error->all(FLERR,"Illegal fix imd parameter: port < 1024"); + if (imd_port < 1024) error->all(FLERR,"Illegal fix imd parameter: port < 1024"); /* default values for optional flags */ - imd_version = 2; + imd_version = 2; unwrap_flag = 0; nowait_flag = 0; connect_msg = 1; imd_fscale = 1.0; imd_trate = 1; - /* IMDv3-only flags. Aren't stored as class attributes since they're converted into IMDSessionInfo */ + // IMDv3-only flags. + // Aren't stored as class attributes since they're converted into IMDSessionInfo + int time_flag = 1; int box_flag = 1; int coord_flag = 1; @@ -507,17 +509,17 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : } else if (0 == strcmp(arg[iarg], "forces")) { force_flag = utils::logical(FLERR, arg[iarg+1], false, lmp); } else { - error->all(FLERR,"Unknown fix imd parameter"); + error->all(FLERR,"Unknown fix imd parameter: {}", arg[iarg]); } ++iarg; ++iarg; } /* sanity check on parameters */ if (imd_trate < 1) - error->all(FLERR,"Illegal fix imd parameter. trate < 1."); + error->all(FLERR,"Illegal fix imd trate parameter. trate < 1."); if (imd_version != 2 && imd_version != 3) - error->all(FLERR,"Illegal fix imd parameter. version != 2 or 3."); + error->all(FLERR,"Illegal fix imd version parameter: version != 2 or 3."); imdsinfo = new IMDSessionInfo; @@ -553,19 +555,11 @@ FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : MPI_Comm_rank(world,&me); /* initialize various imd state variables. */ - clientsock = nullptr; - localsock = nullptr; nlevels_respa = 0; imd_inactive = 0; imd_terminate = 0; imd_forces = 0; - recv_force_buf = nullptr; maxbuf = 0; - coord_data = nullptr; - vel_data = nullptr; - force_data = nullptr; - idmap = nullptr; - rev_idmap = nullptr; if (imd_version == 3) { msglen = 0; @@ -673,8 +667,8 @@ FixIMD::~FixIMD() imdsock_destroy(clientsock); imdsock_shutdown(localsock); imdsock_destroy(localsock); - clientsock=nullptr; - localsock=nullptr; + clientsock = nullptr; + localsock = nullptr; } /* ---------------------------------------------------------------------- */ From 494f58904cb0622fe05527b80e41fb17a5c87121 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 8 Jan 2025 20:54:04 -0500 Subject: [PATCH 13/37] correct reference --- doc/src/fix_imd.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_imd.rst b/doc/src/fix_imd.rst index fd0fd2dc00..a4d7d9d387 100644 --- a/doc/src/fix_imd.rst +++ b/doc/src/fix_imd.rst @@ -71,7 +71,7 @@ that it can monitor the progress of the simulation and interactively apply forces to selected atoms. If LAMMPS is compiled with the pre-processor flag -:ref:`-DLAMMPS_ASYNC_IMD <_misc>` then fix imd will use POSIX threads to +:ref:`-DLAMMPS_ASYNC_IMD ` then fix imd will use POSIX threads to spawn an IMD communication thread on MPI rank 0 in order to offload data exchange with the IMD client from the main execution thread and potentially lower the inferred latencies for slow communication From 2b688bb0131cdbd2b152fd2aec3af25386fc9f45 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 9 Jan 2025 21:32:58 -0500 Subject: [PATCH 14/37] replace defines with static constexpr --- src/QEQ/fix_qeq.cpp | 3 +++ src/QEQ/fix_qeq.h | 7 +------ src/QEQ/fix_qeq_shielded.cpp | 2 ++ 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp index cd331b5f30..bbfb4d8a0b 100644 --- a/src/QEQ/fix_qeq.cpp +++ b/src/QEQ/fix_qeq.cpp @@ -41,6 +41,9 @@ using namespace LAMMPS_NS; using namespace FixConst; static constexpr double QSUMSMALL = 0.00001; +static constexpr int MIN_CAP = 50; +static constexpr double SAFE_ZONE = 1.2; +static constexpr bigint MIN_NBRS = 100; namespace { class qeq_parser_error : public std::exception { diff --git a/src/QEQ/fix_qeq.h b/src/QEQ/fix_qeq.h index 5d0e05c7c5..0124a6ee8a 100644 --- a/src/QEQ/fix_qeq.h +++ b/src/QEQ/fix_qeq.h @@ -16,12 +16,6 @@ #include "fix.h" -#define EV_TO_KCAL_PER_MOL 14.4 -#define DANGER_ZONE 0.90 -#define MIN_CAP 50 -#define SAFE_ZONE 1.2 -#define MIN_NBRS 100 - namespace LAMMPS_NS { class FixQEq : public Fix { @@ -36,6 +30,7 @@ class FixQEq : public Fix { void min_pre_force(int) override; double compute_scalar() override; + static constexpr double DANGER_ZONE = 0.90; // derived child classes must provide these functions diff --git a/src/QEQ/fix_qeq_shielded.cpp b/src/QEQ/fix_qeq_shielded.cpp index b4827d8694..8d1ec7f87e 100644 --- a/src/QEQ/fix_qeq_shielded.cpp +++ b/src/QEQ/fix_qeq_shielded.cpp @@ -33,6 +33,8 @@ using namespace LAMMPS_NS; +static constexpr double EV_TO_KCAL_PER_MOL = 14.4; + /* ---------------------------------------------------------------------- */ FixQEqShielded::FixQEqShielded(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg, arg) From 3731513576ffd97b581ab941d0bb8486182acc86 Mon Sep 17 00:00:00 2001 From: ljwoods2 Date: Thu, 9 Jan 2025 21:55:44 -0700 Subject: [PATCH 15/37] example IMDv3 usage --- examples/PACKAGES/imd/deca-ala-solv_imd_v3.py | 13 ++++++++ examples/PACKAGES/imd/in.deca-ala-solv_imd_v3 | 31 +++++++++++++++++++ 2 files changed, 44 insertions(+) create mode 100644 examples/PACKAGES/imd/deca-ala-solv_imd_v3.py create mode 100644 examples/PACKAGES/imd/in.deca-ala-solv_imd_v3 diff --git a/examples/PACKAGES/imd/deca-ala-solv_imd_v3.py b/examples/PACKAGES/imd/deca-ala-solv_imd_v3.py new file mode 100644 index 0000000000..176e15e811 --- /dev/null +++ b/examples/PACKAGES/imd/deca-ala-solv_imd_v3.py @@ -0,0 +1,13 @@ +""" +For use with 'in.deca-ala-solv_imd_v3'. + +Tested with imdclient v0.1.4 and MDAnalysis v2.8.0 +""" +from imdclient.IMD import IMDReader +import MDAnalysis as mda + +u = mda.Universe('data.deca-ala-solv', "imd://localhost:5678", topology_format='DATA') + +for ts in u.trajectory: + print(ts.time) + print(ts.velocities) \ No newline at end of file diff --git a/examples/PACKAGES/imd/in.deca-ala-solv_imd_v3 b/examples/PACKAGES/imd/in.deca-ala-solv_imd_v3 new file mode 100644 index 0000000000..a7e8244cdc --- /dev/null +++ b/examples/PACKAGES/imd/in.deca-ala-solv_imd_v3 @@ -0,0 +1,31 @@ +# +units real +neighbor 2.5 bin +neigh_modify delay 1 every 1 + +atom_style full +bond_style harmonic +angle_style charmm +dihedral_style charmm +improper_style harmonic + +pair_style lj/charmm/coul/long 8 10 +pair_modify mix arithmetic +special_bonds charmm +read_data data.deca-ala-solv + + +group peptide id <= 103 +fix rigidh all shake 1e-6 100 1000 t 1 2 3 4 5 a 23 + +thermo 100 +thermo_style multi +timestep 2.0 +kspace_style pppm 1e-5 + +fix ensemble all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 0.2 + +# IMD setup. Client code available in 'deca-ala-solv_imd_v3.py' +fix comm all imd 5678 unwrap on trate 10 version 3 time on box on coordinates on velocities on forces off + +run 5000000 \ No newline at end of file From 4f6c3d12f7c59cfde2db6eed6eca498ee5629fcd Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 10 Jan 2025 20:10:53 -0500 Subject: [PATCH 16/37] avoid variable definition error: expected 3 arguments but found N --- src/variable.cpp | 48 +++++++++++++++++++++++++++++++++--------------- 1 file changed, 33 insertions(+), 15 deletions(-) diff --git a/src/variable.cpp b/src/variable.cpp index 135882fde9..031709166b 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -457,15 +457,21 @@ void Variable::set(int narg, char **arg) // data = 2 values, 1st is string to eval, 2nd is filled on retrieval } else if (strcmp(arg[1],"equal") == 0) { - if (narg != 3) - error->all(FLERR,"Illegal variable command: expected 3 arguments but found {}{}", - narg, utils::errorurl(3)); + if (narg < 3) utils::missing_cmd_args(FLERR, "variable equal", error); + + // combine excess arguments into single string with a blank as separator + std::string combined = arg[2]; + for (int iarg = 3; iarg < narg; ++iarg) { + combined += ' '; + combined += arg[iarg]; + } + int ivar = find(arg[0]); if (ivar >= 0) { if (style[ivar] != EQUAL) error->all(FLERR,"Cannot redefine variable as a different style"); delete[] data[ivar][0]; - data[ivar][0] = utils::strdup(arg[2]); + data[ivar][0] = utils::strdup(combined); replaceflag = 1; } else { if (nvar == maxvar) grow(); @@ -474,7 +480,7 @@ void Variable::set(int narg, char **arg) which[nvar] = 0; pad[nvar] = 0; data[nvar] = new char*[num[nvar]]; - data[nvar][0] = utils::strdup(arg[2]); + data[nvar][0] = utils::strdup(combined); data[nvar][1] = new char[VALUELENGTH]; strcpy(data[nvar][1],"(undefined)"); } @@ -485,15 +491,21 @@ void Variable::set(int narg, char **arg) // data = 1 value, string to eval } else if (strcmp(arg[1],"atom") == 0) { - if (narg != 3) - error->all(FLERR,"Illegal variable command: expected 3 arguments but found {}{}", - narg, utils::errorurl(3)); + if (narg < 3) utils::missing_cmd_args(FLERR, "variable atom", error); + + // combine excess arguments into single string with a blank as separator + std::string combined = arg[2]; + for (int iarg = 3; iarg < narg; ++iarg) { + combined += ' '; + combined += arg[iarg]; + } + int ivar = find(arg[0]); if (ivar >= 0) { if (style[ivar] != ATOM) error->all(FLERR,"Cannot redefine variable as a different style"); delete[] data[ivar][0]; - data[ivar][0] = utils::strdup(arg[2]); + data[ivar][0] = utils::strdup(combined); replaceflag = 1; } else { if (nvar == maxvar) grow(); @@ -502,7 +514,7 @@ void Variable::set(int narg, char **arg) which[nvar] = 0; pad[nvar] = 0; data[nvar] = new char*[num[nvar]]; - data[nvar][0] = utils::strdup(arg[2]); + data[nvar][0] = utils::strdup(combined); } // VECTOR @@ -513,16 +525,22 @@ void Variable::set(int narg, char **arg) // immediately store it as N-length vector and set dynamic flag to 0 } else if (strcmp(arg[1],"vector") == 0) { - if (narg != 3) - error->all(FLERR,"Illegal variable command: expected 3 arguments but found {}{}", - narg, utils::errorurl(3)); + if (narg < 3) utils::missing_cmd_args(FLERR, "variable atom", error); + + // combine excess arguments into single string with a blank as separator + std::string combined = arg[2]; + for (int iarg = 3; iarg < narg; ++iarg) { + combined += ' '; + combined += arg[iarg]; + } + int ivar = find(arg[0]); if (ivar >= 0) { if (style[ivar] != VECTOR) error->all(FLERR,"Cannot redefine variable as a different style"); delete[] data[ivar][0]; delete[] data[ivar][1]; - data[ivar][0] = utils::strdup(arg[2]); + data[ivar][0] = utils::strdup(combined); if (data[ivar][0][0] != '[') vecs[ivar].dynamic = 1; else { @@ -539,7 +557,7 @@ void Variable::set(int narg, char **arg) which[nvar] = 0; pad[nvar] = 0; data[nvar] = new char*[num[nvar]]; - data[nvar][0] = utils::strdup(arg[2]); + data[nvar][0] = utils::strdup(combined); if (data[nvar][0][0] != '[') { vecs[nvar].dynamic = 1; data[nvar][1] = nullptr; From c4a59063827f7dc80258073aacaeacce2d321dbb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 10 Jan 2025 20:12:09 -0500 Subject: [PATCH 17/37] stop checking ML-PACE with coverity scan as recent CMake changes break it --- .github/workflows/coverity.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/coverity.yml b/.github/workflows/coverity.yml index c0c3e3f89a..2691b9e895 100644 --- a/.github/workflows/coverity.yml +++ b/.github/workflows/coverity.yml @@ -67,7 +67,6 @@ jobs: -D PKG_MANIFOLD=on \ -D PKG_MDI=on \ -D PKG_MGPT=on \ - -D PKG_ML-PACE=on \ -D PKG_ML-RANN=on \ -D PKG_MOLFILE=on \ -D PKG_NETCDF=on \ From 3126482c4803b898ee98a36a3d3e1c3b52875164 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 10 Jan 2025 20:53:24 -0500 Subject: [PATCH 18/37] update variable command tests to include whitespace --- unittest/commands/test_variables.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/unittest/commands/test_variables.cpp b/unittest/commands/test_variables.cpp index 466dac17cc..0d8d28426c 100644 --- a/unittest/commands/test_variables.cpp +++ b/unittest/commands/test_variables.cpp @@ -140,8 +140,8 @@ TEST_F(VariableTest, CreateDelete) command("variable ten1 universe 1 2 3 4"); command("variable ten2 uloop 4"); command("variable ten3 uloop 4 pad"); - command("variable ten4 vector [0,1,2,3,5,7,11]"); - command("variable ten5 vector [0.5,1.25]"); + command("variable ten4 vector [0,1, 2,3, 5,7,11]"); + command("variable ten5 vector [ 0.5, 1.25 ]"); command("variable dummy index 0"); command("variable file equal is_file(MYFILE)"); command("variable iswin equal is_os(^Windows)"); @@ -323,13 +323,13 @@ TEST_F(VariableTest, Expressions) BEGIN_HIDE_OUTPUT(); command("variable one index 1"); command("variable two equal 2"); - command("variable three equal v_one+v_two"); + command("variable three equal v_one + v_two"); command("variable four equal PI"); command("variable five equal version"); command("variable six equal XXX"); command("variable seven equal -v_one"); command("variable eight equal v_three-0.5"); - command("variable nine equal v_two*(v_one+v_three)"); + command("variable nine equal v_two * (v_one+v_three)"); command("variable ten equal (1.0/v_two)^2"); command("variable eleven equal v_three%2"); command("variable twelve equal 1==2"); @@ -341,7 +341,7 @@ TEST_F(VariableTest, Expressions) command("variable ten8 equal 1|^0"); command("variable ten9 equal v_one-v_ten9"); command("variable ten10 internal 100.0"); - command("variable ten11 equal (1!=1)+(2<1)+(2<=1)+(1>2)+(1>=2)+(1&&0)+(0||0)+(1|^1)+10^0"); + command("variable ten11 equal (1 != 1)+(2 < 1)+(2<=1)+(1>2)+(1>=2)+(1&&0)+(0||0)+(1|^1)+10^0"); command("variable ten12 equal yes+no+on+off+true+false"); command("variable err1 equal v_one/v_ten7"); command("variable err2 equal v_one%v_ten7"); @@ -350,7 +350,7 @@ TEST_F(VariableTest, Expressions) command("variable vec2 vector v_vec1*0.5"); command("variable vec3 equal v_vec2[3]"); command("variable vec4 vector '[1, 5, 2.5, -10, -5, 20, 120, 4, 3, 3]'"); - command("variable sort vector sort(v_vec4)"); + command("variable sort vector sort(v_vec4 )"); command("variable rsrt vector rsort(v_vec4)"); command("variable max2 equal sort(v_vec4)[2]"); command("variable rmax equal rsort(v_vec4)[1]"); From 0d85c5c704afe078350f18a4a2b2f28b24a4ac25 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 11 Jan 2025 11:13:47 -0500 Subject: [PATCH 19/37] add separators --- src/MISC/fix_imd.cpp | 59 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/src/MISC/fix_imd.cpp b/src/MISC/fix_imd.cpp index 493c3a4e0f..0dbafe49cb 100644 --- a/src/MISC/fix_imd.cpp +++ b/src/MISC/fix_imd.cpp @@ -80,6 +80,8 @@ negotiate an appropriate license for such distribution." using namespace LAMMPS_NS; using namespace FixConst; +/* ---------------------------------------------------------------------- */ + /* re-usable integer hash table code with static linkage. */ /** hash table top level data structure */ @@ -672,6 +674,7 @@ FixIMD::~FixIMD() } /* ---------------------------------------------------------------------- */ + int FixIMD::setmask() { int mask = 0; @@ -682,6 +685,7 @@ int FixIMD::setmask() } /* ---------------------------------------------------------------------- */ + void FixIMD::init() { if (utils::strmatch(update->integrate_style,"^respa")) @@ -771,6 +775,8 @@ void FixIMD::setup(int) } } +/* ---------------------------------------------------------------------- */ + void FixIMD::setup_v2() { /* nme: number of atoms in group on this MPI task * nmax: max number of atoms in group across all MPI tasks @@ -858,6 +864,8 @@ void FixIMD::setup_v2() { } +/* ---------------------------------------------------------------------- */ + void FixIMD::setup_v3() { /* nme: number of atoms in group on this MPI task @@ -1026,12 +1034,15 @@ void FixIMD::post_force(int /*vflag*/) } /* ---------------------------------------------------------------------- */ + void FixIMD::post_force_respa(int vflag, int ilevel, int /*iloop*/) { /* only process IMD on the outmost RESPA level. */ if (ilevel == nlevels_respa-1) post_force(vflag); } +/* ---------------------------------------------------------------------- */ + void FixIMD::end_of_step() { if (imd_version == 3 && update->ntimestep % imd_trate == 0) { @@ -1041,11 +1052,14 @@ void FixIMD::end_of_step() /* ---------------------------------------------------------------------- */ /* local memory usage. approximately. */ + double FixIMD::memory_usage() { return static_cast(num_coords+maxbuf+imd_forces)*size_one; } +/* ---------------------------------------------------------------------- */ + void FixIMD::handle_step_v2() { /* check for reconnect */ @@ -1358,6 +1372,8 @@ void FixIMD::handle_step_v2() { } +/* ---------------------------------------------------------------------- */ + void FixIMD::handle_client_input_v3() { // IMDV3 @@ -1542,6 +1558,8 @@ void FixIMD::handle_client_input_v3() { } } +/* ---------------------------------------------------------------------- */ + void FixIMD::handle_output_v3() { tagint *tag = atom->tag; @@ -1818,6 +1836,7 @@ int imdsock_init() { #endif } +/* ---------------------------------------------------------------------- */ void * imdsock_create() { imdsocket * s; @@ -1836,6 +1855,8 @@ void * imdsock_create() { return (void *) s; } +/* ---------------------------------------------------------------------- */ + int imdsock_bind(void * v, int port) { auto s = (imdsocket *) v; auto *addr = &(s->addr); @@ -1847,11 +1868,15 @@ int imdsock_bind(void * v, int port) { return bind(s->sd, (struct sockaddr *) addr, s->addrlen); } +/* ---------------------------------------------------------------------- */ + int imdsock_listen(void * v) { auto s = (imdsocket *) v; return listen(s->sd, 5); } +/* ---------------------------------------------------------------------- */ + void *imdsock_accept(void * v) { int rc; imdsocket *new_s = nullptr, *s = (imdsocket *) v; @@ -1881,6 +1906,8 @@ void *imdsock_accept(void * v) { return (void *)new_s; } +/* ---------------------------------------------------------------------- */ + int imdsock_write(void * v, const void *buf, int len) { auto s = (imdsocket *) v; #if defined(_MSC_VER) || defined(__MINGW32__) @@ -1890,6 +1917,8 @@ int imdsock_write(void * v, const void *buf, int len) { #endif } +/* ---------------------------------------------------------------------- */ + int imdsock_read(void * v, void *buf, int len) { auto s = (imdsocket *) v; #if defined(_MSC_VER) || defined(__MINGW32__) @@ -1900,6 +1929,8 @@ int imdsock_read(void * v, void *buf, int len) { } +/* ---------------------------------------------------------------------- */ + void imdsock_shutdown(void *v) { auto s = (imdsocket *) v; if (s == nullptr) @@ -1912,6 +1943,8 @@ void imdsock_shutdown(void *v) { #endif } +/* ---------------------------------------------------------------------- */ + void imdsock_destroy(void * v) { auto s = (imdsocket *) v; if (s == nullptr) @@ -1925,6 +1958,8 @@ void imdsock_destroy(void * v) { free(s); } +/* ---------------------------------------------------------------------- */ + int imdsock_selread(void *v, int sec) { auto s = (imdsocket *)v; fd_set rfd; @@ -1944,6 +1979,8 @@ int imdsock_selread(void *v, int sec) { } +/* ---------------------------------------------------------------------- */ + int imdsock_selwrite(void *v, int sec) { auto s = (imdsocket *)v; fd_set wfd; @@ -1979,6 +2016,8 @@ typedef union { } b; } netint; +/* ---------------------------------------------------------------------- */ + static int32 imd_htonl(int32 h) { netint n; n.b.highest = h >> 24; @@ -1988,22 +2027,30 @@ static int32 imd_htonl(int32 h) { return n.i; } +/* ---------------------------------------------------------------------- */ + static int32 imd_ntohl(int32 n) { netint u; u.i = n; return (u.b.highest << 24 | u.b.high << 16 | u.b.low << 8 | u.b.lowest); } +/* ---------------------------------------------------------------------- */ + static void imd_fill_header(IMDheader *header, IMDType type, int32 length) { header->type = imd_htonl((int32)type); header->length = imd_htonl(length); } +/* ---------------------------------------------------------------------- */ + static void swap_header(IMDheader *header) { header->type = imd_ntohl(header->type); header->length= imd_ntohl(header->length); } +/* ---------------------------------------------------------------------- */ + static int32 imd_readn(void *s, char *ptr, int32 n) { int32 nleft; int32 nread; @@ -2023,6 +2070,8 @@ static int32 imd_readn(void *s, char *ptr, int32 n) { return n-nleft; } +/* ---------------------------------------------------------------------- */ + static int32 imd_writen(void *s, const char *ptr, int32 n) { int32 nleft; int32 nwritten; @@ -2041,6 +2090,8 @@ static int32 imd_writen(void *s, const char *ptr, int32 n) { return n; } +/* ---------------------------------------------------------------------- */ + int imd_handshake_v2(void *s) { IMDheader header; imd_fill_header(&header, IMD_HANDSHAKE, 1); @@ -2048,6 +2099,8 @@ int imd_handshake_v2(void *s) { return (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE); } +/* ---------------------------------------------------------------------- */ + int imd_handshake_v3(void *s, IMDSessionInfo *imdsinfo) { IMDheader header; imd_fill_header(&header, IMD_HANDSHAKE, 1); @@ -2081,17 +2134,23 @@ IMDType imd_recv_header(void *s, int32 *length) { return IMDType(header.type); } +/* ---------------------------------------------------------------------- */ + int imd_recv_mdcomm(void *s, int32 n, int32 *indices, float *forces) { if (imd_readn(s, (char *)indices, 4*n) != 4*n) return 1; if (imd_readn(s, (char *)forces, 12*n) != 12*n) return 1; return 0; } +/* ---------------------------------------------------------------------- */ + int imd_recv_energies(void *s, IMDEnergies *energies) { return (imd_readn(s, (char *)energies, sizeof(IMDEnergies)) != sizeof(IMDEnergies)); } +/* ---------------------------------------------------------------------- */ + int imd_recv_fcoords(void *s, int32 n, float *coords) { return (imd_readn(s, (char *)coords, 12*n) != 12*n); } From b508229bd2259a11c82ff2e077b8d5da29029a29 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 11 Jan 2025 18:05:16 -0500 Subject: [PATCH 20/37] adapt to LAMMPS programming style --- src/fix_deposit.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp index d6630bcbf0..c9d09fa8b7 100644 --- a/src/fix_deposit.cpp +++ b/src/fix_deposit.cpp @@ -208,14 +208,14 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) : FixDeposit::~FixDeposit() { delete random; - delete [] molfrac; - delete [] idrigid; - delete [] idshake; - delete [] idregion; - delete [] vstr; - delete [] xstr; - delete [] ystr; - delete [] zstr; + delete[] molfrac; + delete[] idrigid; + delete[] idshake; + delete[] idregion; + delete[] vstr; + delete[] xstr; + delete[] ystr; + delete[] zstr; memory->destroy(coords); memory->destroy(imageflags); } @@ -737,7 +737,7 @@ void FixDeposit::options(int narg, char **arg) mode = MOLECULE; onemols = &atom->molecules[imol]; nmol = onemols[0]->nset; - delete [] molfrac; + delete[] molfrac; molfrac = new double[nmol]; molfrac[0] = 1.0/nmol; for (int i = 1; i < nmol-1; i++) molfrac[i] = molfrac[i-1] + 1.0/nmol; @@ -755,13 +755,13 @@ void FixDeposit::options(int narg, char **arg) iarg += nmol+1; } else if (strcmp(arg[iarg],"rigid") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - delete [] idrigid; + delete[] idrigid; idrigid = utils::strdup(arg[iarg+1]); rigidflag = 1; iarg += 2; } else if (strcmp(arg[iarg],"shake") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command"); - delete [] idshake; + delete[] idshake; idshake = utils::strdup(arg[iarg+1]); shakeflag = 1; iarg += 2; From 61c541ff929a7d4bfb7dadcf89d1f9711f434ea1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 11 Jan 2025 18:47:35 -0500 Subject: [PATCH 21/37] add github action to build LAMMPS-GUI as flatpak --- .github/workflows/lammps-gui-flatpak.yml | 52 ++++++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 .github/workflows/lammps-gui-flatpak.yml diff --git a/.github/workflows/lammps-gui-flatpak.yml b/.github/workflows/lammps-gui-flatpak.yml new file mode 100644 index 0000000000..7fe6b3180b --- /dev/null +++ b/.github/workflows/lammps-gui-flatpak.yml @@ -0,0 +1,52 @@ +# GitHub action to build LAMMPS-GUI as a flatpak bundle +name: "Build LAMMPS-GUI as flatpak bundle" + +on: + push: + branches: + - develop + - collected-small-changes + pull_request: + branches: + - develop + + workflow_dispatch: + +concurrency: + group: ${{ github.event_name }}-${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{github.event_name == 'pull_request'}} + +jobs: + build: + name: Linux flatpak build + if: ${{ github.repository == 'akohlmey/lammps' }} + runs-on: ubuntu-latest + + steps: + - name: Checkout repository + uses: actions/checkout@v4 + with: + fetch-depth: 2 + + - name: Install extra packages + run: | + sudo apt-get update + sudo apt-get install -y ccache \ + libeigen3-dev \ + libcurl4-openssl-dev \ + mold \ + ninja-build \ + python3-dev \ + flatpak \ + flatpak-builder + + - name: Set up access to flatpak repo + run: flatpak --user remote-add --if-not-exists flathub https://dl.flathub.org/repo/flathub.flatpakrepo + + - name: Build flatpak + run: | + mkdir flatpack-state + flatpak-builder --force-clean --verbose --repo=flatpak-repo + --install-deps-from=flathub --state-dir=flatpak-state + --user --ccache --default-branch=${{ github.ref_name }} + flatpak-build tools/lammps-gui/org.lammps.lammps-gui.yml From cb2acb633b75d9544fbcbc5a505530c3f1d77539 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 11 Jan 2025 18:56:59 -0500 Subject: [PATCH 22/37] update workflow --- .github/workflows/lammps-gui-flatpak.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/lammps-gui-flatpak.yml b/.github/workflows/lammps-gui-flatpak.yml index 7fe6b3180b..2e3cd7eb2f 100644 --- a/.github/workflows/lammps-gui-flatpak.yml +++ b/.github/workflows/lammps-gui-flatpak.yml @@ -46,7 +46,7 @@ jobs: - name: Build flatpak run: | mkdir flatpack-state - flatpak-builder --force-clean --verbose --repo=flatpak-repo - --install-deps-from=flathub --state-dir=flatpak-state - --user --ccache --default-branch=${{ github.ref_name }} + flatpak-builder --force-clean --verbose --repo=flatpak-repo \ + --install-deps-from=flathub --state-dir=flatpak-state \ + --user --ccache --default-branch=${{ github.ref_name }} \ flatpak-build tools/lammps-gui/org.lammps.lammps-gui.yml From fb1c090f19aba2d9a4aa6ac9ead44ac7dc68fa52 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 11 Jan 2025 18:59:22 -0500 Subject: [PATCH 23/37] revert to lammps repo --- .github/workflows/lammps-gui-flatpak.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/lammps-gui-flatpak.yml b/.github/workflows/lammps-gui-flatpak.yml index 2e3cd7eb2f..d7a2965a04 100644 --- a/.github/workflows/lammps-gui-flatpak.yml +++ b/.github/workflows/lammps-gui-flatpak.yml @@ -19,7 +19,7 @@ concurrency: jobs: build: name: Linux flatpak build - if: ${{ github.repository == 'akohlmey/lammps' }} + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest steps: From 672faf9fe5b20df4e48619422c798277d055b3fa Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 11 Jan 2025 19:06:49 -0500 Subject: [PATCH 24/37] operate on "develop" branch --- .github/workflows/lammps-gui-flatpak.yml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/.github/workflows/lammps-gui-flatpak.yml b/.github/workflows/lammps-gui-flatpak.yml index d7a2965a04..b5e1f169d8 100644 --- a/.github/workflows/lammps-gui-flatpak.yml +++ b/.github/workflows/lammps-gui-flatpak.yml @@ -6,9 +6,6 @@ on: branches: - develop - collected-small-changes - pull_request: - branches: - - develop workflow_dispatch: @@ -46,6 +43,7 @@ jobs: - name: Build flatpak run: | mkdir flatpack-state + sed -i -e 's/branch:.*/branch: develop/' tools/lammps-gui/org.lammps.lammps-gui.yml flatpak-builder --force-clean --verbose --repo=flatpak-repo \ --install-deps-from=flathub --state-dir=flatpak-state \ --user --ccache --default-branch=${{ github.ref_name }} \ From cf2e800aaaeaee58eda5e42b402d9ecfe7f4a799 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 11 Jan 2025 19:20:06 -0500 Subject: [PATCH 25/37] extract bundle from local flatpak repo and try to install it --- .github/workflows/lammps-gui-flatpak.yml | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/.github/workflows/lammps-gui-flatpak.yml b/.github/workflows/lammps-gui-flatpak.yml index b5e1f169d8..6e2077b229 100644 --- a/.github/workflows/lammps-gui-flatpak.yml +++ b/.github/workflows/lammps-gui-flatpak.yml @@ -16,7 +16,7 @@ concurrency: jobs: build: name: Linux flatpak build - if: ${{ github.repository == 'lammps/lammps' }} + if: ${{ github.repository == 'akohlmey/lammps' }} runs-on: ubuntu-latest steps: @@ -48,3 +48,7 @@ jobs: --install-deps-from=flathub --state-dir=flatpak-state \ --user --ccache --default-branch=${{ github.ref_name }} \ flatpak-build tools/lammps-gui/org.lammps.lammps-gui.yml + flatpak build-bundle --runtime-repo=https://flathub.org/repo/flathub.flatpakrepo \ + --verbose flatpak-repo LAMMPS-Linux-x86_64-GUI.flatpak \ + org.lammps.lammps-gui ${{ github.ref_name }} + flatpak install -y -v --user LAMMPS-Linux-x86_64-GUI.flatpak From be596cca0c2a4f661c55b434dca0f93ace8790d8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 11 Jan 2025 19:31:31 -0500 Subject: [PATCH 26/37] revert setting to run flatpak builder test action only after merges --- .github/workflows/lammps-gui-flatpak.yml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/lammps-gui-flatpak.yml b/.github/workflows/lammps-gui-flatpak.yml index 6e2077b229..d7dc602476 100644 --- a/.github/workflows/lammps-gui-flatpak.yml +++ b/.github/workflows/lammps-gui-flatpak.yml @@ -5,7 +5,6 @@ on: push: branches: - develop - - collected-small-changes workflow_dispatch: @@ -15,8 +14,8 @@ concurrency: jobs: build: - name: Linux flatpak build - if: ${{ github.repository == 'akohlmey/lammps' }} + name: LAMMPS-GUI flatpak build + if: ${{ github.repository == 'lammps/lammps' }} runs-on: ubuntu-latest steps: From dce1231052f7329a4fa211d51d5ca655266dbd9c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 12 Jan 2025 00:39:18 -0500 Subject: [PATCH 27/37] add a lammps_eval() function to the C library interface and all derived wrappers --- doc/src/Fortran.rst | 2 ++ doc/src/Library_objects.rst | 6 +++++ examples/COUPLE/plugin/liblammpsplugin.c | 1 + examples/COUPLE/plugin/liblammpsplugin.h | 1 + fortran/lammps.f90 | 24 +++++++++++++++++-- python/lammps/core.py | 27 +++++++++++++++++++++ src/library.cpp | 30 ++++++++++++++++++++++++ src/library.h | 1 + tools/swig/lammps.i | 4 ++++ 9 files changed, 94 insertions(+), 2 deletions(-) diff --git a/doc/src/Fortran.rst b/doc/src/Fortran.rst index d50c7bdcf3..bc641a237f 100644 --- a/doc/src/Fortran.rst +++ b/doc/src/Fortran.rst @@ -321,6 +321,8 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS. :ftype set_string_variable: subroutine :f set_internal_variable: :f:subr:`set_internal_variable` :ftype set_internal_variable: subroutine + :f eval: :f:func:`eval` + :ftype eval: function :f gather_atoms: :f:subr:`gather_atoms` :ftype gather_atoms: subroutine :f gather_atoms_concat: :f:subr:`gather_atoms_concat` diff --git a/doc/src/Library_objects.rst b/doc/src/Library_objects.rst index 7c0ca824d7..89056a5db0 100644 --- a/doc/src/Library_objects.rst +++ b/doc/src/Library_objects.rst @@ -12,6 +12,7 @@ fixes, or variables in LAMMPS using the following functions: - :cpp:func:`lammps_set_string_variable` - :cpp:func:`lammps_set_internal_variable` - :cpp:func:`lammps_variable_info` +- :cpp:func:`lammps_eval` ----------------------- @@ -55,6 +56,11 @@ fixes, or variables in LAMMPS using the following functions: ----------------------- +.. doxygenfunction:: lammps_eval + :project: progguide + +----------------------- + .. doxygenenum:: _LMP_DATATYPE_CONST .. doxygenenum:: _LMP_STYLE_CONST diff --git a/examples/COUPLE/plugin/liblammpsplugin.c b/examples/COUPLE/plugin/liblammpsplugin.c index 5b1308a5cc..99e38df32b 100644 --- a/examples/COUPLE/plugin/liblammpsplugin.c +++ b/examples/COUPLE/plugin/liblammpsplugin.c @@ -117,6 +117,7 @@ liblammpsplugin_t *liblammpsplugin_load(const char *lib) ADDSYM(set_string_variable); ADDSYM(set_internal_variable); ADDSYM(variable_info); + ADDSYM(eval); ADDSYM(gather_atoms); ADDSYM(gather_atoms_concat); diff --git a/examples/COUPLE/plugin/liblammpsplugin.h b/examples/COUPLE/plugin/liblammpsplugin.h index b2df98d630..dd1c9628f5 100644 --- a/examples/COUPLE/plugin/liblammpsplugin.h +++ b/examples/COUPLE/plugin/liblammpsplugin.h @@ -163,6 +163,7 @@ struct _liblammpsplugin { int (*set_string_variable)(void *, const char *, const char *); int (*set_internal_variable)(void *, const char *, double); int (*variable_info)(void *, int, char *, int); + double (*eval)(void *, const char *); void (*gather_atoms)(void *, const char *, int, int, void *); void (*gather_atoms_concat)(void *, const char *, int, int, void *); diff --git a/fortran/lammps.f90 b/fortran/lammps.f90 index 1cc6992c27..552b3dfad3 100644 --- a/fortran/lammps.f90 +++ b/fortran/lammps.f90 @@ -126,6 +126,7 @@ MODULE LIBLAMMPS PROCEDURE :: set_variable => lmp_set_variable PROCEDURE :: set_string_variable => lmp_set_string_variable PROCEDURE :: set_internal_variable => lmp_set_internal_variable + PROCEDURE :: eval => lmp_eval PROCEDURE, PRIVATE :: lmp_gather_atoms_int PROCEDURE, PRIVATE :: lmp_gather_atoms_double GENERIC :: gather_atoms => lmp_gather_atoms_int, & @@ -618,7 +619,14 @@ MODULE LIBLAMMPS INTEGER(c_int) :: lammps_set_internal_variable END FUNCTION lammps_set_internal_variable - SUBROUTINE lammps_gather_atoms(handle, name, type, count, data) BIND(C) + FUNCTION lammps_eval(handle, expr) BIND(C) + IMPORT :: c_ptr, c_double + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, expr + REAL(c_double) :: lammps_eval + END FUNCTION lammps_eval + + SUBROUTINE lammps_gather_atoms(handle, name, TYPE, count, DATA) BIND(C) IMPORT :: c_int, c_ptr IMPLICIT NONE TYPE(c_ptr), VALUE :: handle, name, data @@ -1812,7 +1820,7 @@ CONTAINS SUBROUTINE lmp_set_internal_variable(self, name, val) CLASS(lammps), INTENT(IN) :: self CHARACTER(LEN=*), INTENT(IN) :: name - REAL(KIND=c_double), INTENT(IN) :: val + REAL(c_double), INTENT(IN) :: val INTEGER :: err TYPE(c_ptr) :: Cname @@ -1826,6 +1834,18 @@ CONTAINS END IF END SUBROUTINE lmp_set_internal_variable + ! equivalent function to lammps_eval + FUNCTION lmp_eval(self, expr) + CLASS(lammps), INTENT(IN) :: self + CHARACTER(LEN=*), INTENT(IN) :: expr + REAL(c_double) :: lmp_eval + TYPE(c_ptr) :: Cexpr + + Cexpr = f2c_string(expr) + lmp_eval = lammps_eval(self%handle, Cexpr) + CALL lammps_free(Cexpr) + END FUNCTION lmp_eval + ! equivalent function to lammps_gather_atoms (for integers) SUBROUTINE lmp_gather_atoms_int(self, name, count, data) CLASS(lammps), INTENT(IN) :: self diff --git a/python/lammps/core.py b/python/lammps/core.py index 30df44f050..d6fb63b815 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -339,6 +339,9 @@ class lammps(object): self.lib.lammps_extract_variable_datatype.argtypes = [c_void_p, c_char_p] self.lib.lammps_extract_variable_datatype.restype = c_int + self.lib.lammps_eval.argtypes = [c_void_p, c_char_p] + self.lib.lammps_eval.restype = c_double + self.lib.lammps_fix_external_get_force.argtypes = [c_void_p, c_char_p] self.lib.lammps_fix_external_get_force.restype = POINTER(POINTER(c_double)) @@ -1533,6 +1536,30 @@ class lammps(object): # ------------------------------------------------------------------------- + def eval(self, expr): + """ Evaluate a LAMMPS immediate variable expression + + .. versionadded:: TBD + + This function is a wrapper around the function :cpp:func:`lammps_eval` + of the C library interface. It evaluates and expression like in + immediate variables and returns the value. + + :param expr: immediate variable expression + :type name: string + :return: the result of the evaluation + :rtype: c_double + """ + + if expr: newexpr = expr.encode() + else: return None + + with ExceptionCheck(self): + return self.lib.lammps_eval(self.lmp, newexpr) + return None + + # ------------------------------------------------------------------------- + # return vector of atom properties gathered across procs # 3 variants to match src/library.cpp # name = atom property recognized by LAMMPS in atom->extract() diff --git a/src/library.cpp b/src/library.cpp index cee1304db6..cb5c62aa87 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -2904,6 +2904,36 @@ int lammps_variable_info(void *handle, int idx, char *buffer, int buf_size) { return 0; } +/** Evaluate an immediate variable expression + * +\verbatim embed:rst + +.. versionadded:: TBD + +This function takes a string with an expression, like what can be used +for :doc:`equal style variables `, evaluates it and returns +the resulting (scalar) value as a floating point number. + +\endverbatim + + * \param handle pointer to a previously created LAMMPS instance cast to ``void *``. + * \param expr string with expression + * \return result from expression */ + +double lammps_eval(void *handle, const char *expr) +{ + auto lmp = (LAMMPS *) handle; + double result = 0.0; + + BEGIN_CAPTURE + { + result = lmp->input->variable->compute_equal(expr); + } + END_CAPTURE + + return result; +} + // ---------------------------------------------------------------------- // Library functions for scatter/gather operations of data // ---------------------------------------------------------------------- diff --git a/src/library.h b/src/library.h index 36e67470ae..f64ae9e7c2 100644 --- a/src/library.h +++ b/src/library.h @@ -189,6 +189,7 @@ int lammps_set_variable(void *handle, const char *name, const char *str); int lammps_set_string_variable(void *handle, const char *name, const char *str); int lammps_set_internal_variable(void *handle, const char *name, double value); int lammps_variable_info(void *handle, int idx, char *buf, int bufsize); +double lammps_eval(void *handle, const char *expr); /* ---------------------------------------------------------------------- * Library functions for scatter/gather operations of data diff --git a/tools/swig/lammps.i b/tools/swig/lammps.i index 119edc4367..bc9f46f407 100644 --- a/tools/swig/lammps.i +++ b/tools/swig/lammps.i @@ -141,6 +141,8 @@ extern int lammps_extract_variable_datatype(void *handle, const char *name); extern int lammps_set_variable(void *, const char *, const char *); extern int lammps_set_string_variable(void *, const char *, const char *); extern int lammps_set_internal_variable(void *, const char *, double); +extern int lammps_variable_info(void *handle, int idx, char *buf, int bufsize); +extern double lammps_eval(void *handle, const char *expr); extern void lammps_gather_atoms(void *, char *, int, int, void *); extern void lammps_gather_atoms_concat(void *, char *, int, int, void *); @@ -332,6 +334,8 @@ extern int lammps_extract_variable_datatype(void *handle, const char *name); extern int lammps_set_variable(void *, const char *, const char *); extern int lammps_set_string_variable(void *, const char *, const char *); extern int lammps_set_internal_variable(void *, const char *, double); +extern int lammps_variable_info(void *handle, int idx, char *buf, int bufsize); +extern double lammps_eval(void *handle, const char *expr); extern void lammps_gather_atoms(void *, char *, int, int, void *); extern void lammps_gather_atoms_concat(void *, char *, int, int, void *); From 85dec585a6ff0f1fac0a5e9594d191d27cd67a83 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 12 Jan 2025 14:49:07 -0500 Subject: [PATCH 28/37] add tests for lammps_eval() and its python counterpart --- src/library.cpp | 8 +++++- unittest/c-library/test_library_objects.cpp | 31 +++++++++++++++++++++ unittest/python/python-commands.py | 21 ++++++++++++++ 3 files changed, 59 insertions(+), 1 deletion(-) diff --git a/src/library.cpp b/src/library.cpp index cb5c62aa87..b53257fc45 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -516,6 +516,9 @@ treated as part of the command and will **not** start a second command. The function returns the expanded string in a new string buffer that must be freed with :cpp:func:`lammps_free` after use to avoid a memory leak. +*See also* + :cpp:func:`lammps_eval' + \endverbatim * * \param handle pointer to a previously created LAMMPS instance @@ -2910,10 +2913,13 @@ int lammps_variable_info(void *handle, int idx, char *buffer, int buf_size) { .. versionadded:: TBD -This function takes a string with an expression, like what can be used +This function takes a string with an expression that can be used for :doc:`equal style variables `, evaluates it and returns the resulting (scalar) value as a floating point number. +*See also* + :cpp:func:`lammps_expand' + \endverbatim * \param handle pointer to a previously created LAMMPS instance cast to ``void *``. diff --git a/unittest/c-library/test_library_objects.cpp b/unittest/c-library/test_library_objects.cpp index a175952dfd..4a51381e64 100644 --- a/unittest/c-library/test_library_objects.cpp +++ b/unittest/c-library/test_library_objects.cpp @@ -233,3 +233,34 @@ TEST_F(LibraryObjects, expand) EXPECT_THAT((char *)ptr, StrEq("'xx_$(4+5)_$(PI) ${one}-${two}-${three}'")); lammps_free(ptr); } + +TEST_F(LibraryObjects, eval) +{ + lammps_get_last_error_message(lmp, nullptr, 1); + ::testing::internal::CaptureStdout(); + lammps_commands_string(lmp, "region box1 block 0 10 0 5 -0.5 0.5\n" + "lattice fcc 0.8\n" + "create_box 1 box1\n" + "create_atoms 1 box\n" + "mass * 1.0\n" + "pair_style lj/cut 4.0\n" + "pair_coeff * * 1.0 1.0\n" + "variable t equal 15.0\n" + "velocity all create 1.5 532656\n" + "fix 1 all nve\n" + "run 0 post no\n"); + lammps_command(lmp, "variable one index 1 2 3 4"); + lammps_command(lmp, "variable two equal 2"); + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + ASSERT_EQ(lammps_has_error(lmp), 0); + + EXPECT_DOUBLE_EQ(lammps_eval(lmp, "4+5"), 9.0); + EXPECT_EQ(lammps_has_error(lmp), 0); + EXPECT_DOUBLE_EQ(lammps_eval(lmp, "v_one / 2.0"), 0.5); + EXPECT_EQ(lammps_has_error(lmp), 0); + EXPECT_DOUBLE_EQ(lammps_eval(lmp, "count(all)"), 36.0); + EXPECT_EQ(lammps_has_error(lmp), 0); + EXPECT_DOUBLE_EQ(lammps_eval(lmp, "pe"), -3.9848867644689534); + EXPECT_EQ(lammps_has_error(lmp), 0); +} diff --git a/unittest/python/python-commands.py b/unittest/python/python-commands.py index 2066514318..06dcacd412 100644 --- a/unittest/python/python-commands.py +++ b/unittest/python/python-commands.py @@ -538,6 +538,27 @@ create_atoms 1 single & expanded = self.lmp.expand("'xx_$(4+5)_$(PI) ${one}-${two}-${three}'") self.assertEqual(expanded, "'xx_$(4+5)_$(PI) ${one}-${two}-${three}'") + def test_eval(self): + self.lmp.commands_string( + """region box1 block 0 10 0 5 -0.5 0.5 + lattice fcc 0.8 + create_box 1 box1 + create_atoms 1 box + mass * 1.0 + pair_style lj/cut 4.0 + pair_coeff * * 1.0 1.0 + variable t equal 15.0 + velocity all create 1.5 532656 + fix 1 all nve + run 0 post no""") + self.lmp.command("variable one index 1 2 3 4") + self.lmp.command("variable two equal 2") + + self.assertEqual(self.lmp.eval("4+5"), 9.0) + self.assertEqual(self.lmp.eval("v_one / 2.0"), 0.5) + self.assertEqual(self.lmp.eval("count(all)"), 36.0) + self.assertEqual(self.lmp.eval("pe"), -3.9848867644689534) + def test_get_thermo(self): self.lmp.command("units lj") self.lmp.command("atom_style atomic") From a5c3305c423d1ecdd71d327ece17de9d31d56e5a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 12 Jan 2025 15:09:42 -0500 Subject: [PATCH 29/37] sync with lammps tutorials paper description --- doc/src/Howto_lammps_gui.rst | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/doc/src/Howto_lammps_gui.rst b/doc/src/Howto_lammps_gui.rst index 08be9fb00c..63cf859c57 100644 --- a/doc/src/Howto_lammps_gui.rst +++ b/doc/src/Howto_lammps_gui.rst @@ -64,13 +64,18 @@ simple LAMMPS simulations. It is very suitable for tutorials on LAMMPS since you only need to learn how to use a single program for most tasks and thus time can be saved and people can focus on learning LAMMPS. The tutorials at https://lammpstutorials.github.io/ are specifically -updated for use with LAMMPS-GUI. +updated for use with LAMMPS-GUI and can their tutorial materials can +be downloaded and loaded directly from the GUI. Another design goal is to keep the barrier low when replacing part of the functionality of LAMMPS-GUI with external tools. That said, LAMMPS-GUI has some unique functionality that is not found elsewhere: - auto-adapting to features available in the integrated LAMMPS library +- auto-completion for LAMMPS commands and options +- context-sensitive online help +- start and stop of simulations via mouse or keyboard +- monitoring of simulation progress - interactive visualization using the :doc:`dump image ` command with the option to copy-paste the resulting settings - automatic slide show generation from dump image out at runtime From 0aadc4cf4667b4afbde3b5747533df4413333817 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 12 Jan 2025 23:24:48 -0500 Subject: [PATCH 30/37] use mutex to avoid race condition when accessing thermo data during run --- src/finish.cpp | 2 +- src/library.cpp | 22 ++++++++++-- src/thermo.cpp | 45 ++++++++++++++++++++++--- src/thermo.h | 9 +++++ tools/lammps-gui/lammps-gui.appdata.xml | 2 ++ tools/lammps-gui/lammpsgui.cpp | 2 ++ 6 files changed, 75 insertions(+), 7 deletions(-) diff --git a/src/finish.cpp b/src/finish.cpp index c774f51e58..3b824f2ac4 100644 --- a/src/finish.cpp +++ b/src/finish.cpp @@ -121,9 +121,9 @@ void Finish::end(int flag) MPI_Allreduce(&cpu_loop,&tmp,1,MPI_DOUBLE,MPI_SUM,world); cpu_loop = tmp/nprocs; if (time_loop > 0.0) cpu_loop = cpu_loop/time_loop*100.0; + output->thermo->footer(); if (me == 0) { - output->thermo->footer(); int ntasks = nprocs * nthreads; utils::logmesg(lmp,"Loop time of {:.6g} on {} procs for {} steps with {} atoms\n\n", time_loop,ntasks,update->nsteps,atom->natoms); diff --git a/src/library.cpp b/src/library.cpp index b53257fc45..aa835d0906 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -853,7 +853,11 @@ This differs from :cpp:func:`lammps_get_thermo` in that it does **not** trigger an evaluation. Instead it provides direct access to a read-only location of the last thermo output data and the corresponding keyword strings. How to handle the return value depends on the value of the *what* -argument string. +argument string. When accessing the data from a concurrent thread while +LAMMPS is running, the cache needs to be locked first and then unlocked +after the data is obtained, so that the data is not corrupted while +reading in case LAMMPS wants to update it at the same time. Outside +of a run, the lock/unlock calls have no effect. .. note:: @@ -902,6 +906,14 @@ argument string. - actual field data for column - pointer to int, int64_t or double - yes + * - lock + - acquires lock to thermo data cache + - NULL pointer + - no + * - unlock + - releases lock to thermo data cache + - NULL pointer + - no \endverbatim * @@ -957,8 +969,14 @@ void *lammps_last_thermo(void *handle, const char *what, int index) } else if (field.type == multitype::LAMMPS_DOUBLE) { val = (void *) &field.data.d; } - + } else if (strcmp(what, "lock") == 0) { + th->lock_cache(); + val = nullptr; + } else if (strcmp(what, "unlock") == 0) { + th->unlock_cache(); + val = nullptr; } else val = nullptr; + } END_CAPTURE return val; diff --git a/src/thermo.cpp b/src/thermo.cpp index b2e4f7f934..d7a00fd14c 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -48,6 +48,7 @@ #include #include +#include #include using namespace LAMMPS_NS; @@ -100,8 +101,8 @@ static char fmtbuf[512]; /* ---------------------------------------------------------------------- */ Thermo::Thermo(LAMMPS *_lmp, int narg, char **arg) : - Pointers(_lmp), style(nullptr), vtype(nullptr), field2index(nullptr), argindex1(nullptr), - argindex2(nullptr), temperature(nullptr), pressure(nullptr), pe(nullptr) + Pointers(_lmp), style(nullptr), vtype(nullptr), cache_mutex(nullptr), field2index(nullptr), + argindex1(nullptr), argindex2(nullptr), temperature(nullptr), pressure(nullptr), pe(nullptr) { style = utils::strdup(arg[0]); @@ -208,6 +209,7 @@ void Thermo::init() ValueTokenizer *format_line = nullptr; if (format_line_user.size()) format_line = new ValueTokenizer(format_line_user); + lock_cache(); field_data.clear(); field_data.resize(nfield); std::string format_this, format_line_user_def; @@ -277,6 +279,7 @@ void Thermo::init() format[i] += fmt::format("{:<8} = {} ", keyword[i], format_this); } } + unlock_cache(); // chop off trailing blank or add closing bracket if needed and then add newline if (lineflag == ONELINE) @@ -320,6 +323,9 @@ void Thermo::init() if (index_press_scalar >= 0) pressure = computes[index_press_scalar]; if (index_press_vector >= 0) pressure = computes[index_press_vector]; if (index_pe >= 0) pe = computes[index_pe]; + + // create mutex to protect access to cached thermo data + cache_mutex = new std::mutex; } /* ---------------------------------------------------------------------- @@ -366,9 +372,17 @@ void Thermo::header() /* ---------------------------------------------------------------------- */ +// called at the end of a run from Finish class + void Thermo::footer() { - if (lineflag == YAMLLINE) utils::logmesg(lmp, "...\n"); + if (comm->me == 0) { + if (lineflag == YAMLLINE) utils::logmesg(lmp, "...\n"); + } + + // no more locking for cached thermo data access needed + delete cache_mutex; + cache_mutex = nullptr; } /* ---------------------------------------------------------------------- */ @@ -422,6 +436,7 @@ void Thermo::compute(int flag) } // add each thermo value to line with its specific format + lock_cache(); field_data.clear(); field_data.resize(nfield); @@ -441,6 +456,7 @@ void Thermo::compute(int flag) field_data[ifield] = bivalue; } } + unlock_cache(); // print line to screen and logfile @@ -579,7 +595,8 @@ void Thermo::modify_params(int narg, char **arg) if (iarg + 2 > narg) error->all(FLERR, "Illegal thermo_modify command"); triclinic_general = utils::logical(FLERR, arg[iarg + 1], false, lmp); if (triclinic_general && !domain->triclinic_general) - error->all(FLERR,"Thermo_modify triclinic/general cannot be used " + error->all(FLERR, + "Thermo_modify triclinic/general cannot be used " "if simulation box is not general triclinic"); iarg += 2; @@ -1566,6 +1583,26 @@ int Thermo::evaluate_keyword(const std::string &word, double *answer) return 0; } +/* ---------------------------------------------------------------------- */ + +// lock cache for current thermo data + +void Thermo::lock_cache() +{ + // no locking outside of a run + if (!cache_mutex) return; + cache_mutex->lock(); +} + +// unlock cache for current thermo data + +void Thermo::unlock_cache() +{ + // no locking outside of a run + if (!cache_mutex) return; + cache_mutex->unlock(); +} + /* ---------------------------------------------------------------------- extraction of Compute, Fix, Variable results compute/fix are normalized by atoms if returning extensive value diff --git a/src/thermo.h b/src/thermo.h index 52f69e7609..b5934f48b1 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -17,6 +17,10 @@ #include "pointers.h" #include +namespace std { +class mutex; +} + namespace LAMMPS_NS { class Thermo : protected Pointers { @@ -43,6 +47,8 @@ class Thermo : protected Pointers { int evaluate_keyword(const std::string &, double *); // for accessing cached thermo and related data + void lock_cache(); + void unlock_cache(); const int *get_line() const { return &nline; } const char *get_image_fname() const { return image_fname.c_str(); } @@ -82,6 +88,9 @@ class Thermo : protected Pointers { int nline; std::string image_fname; + // mutex for locking the cache + std::mutex *cache_mutex; + // data used by routines that compute single values int ivalue; // integer value to print diff --git a/tools/lammps-gui/lammps-gui.appdata.xml b/tools/lammps-gui/lammps-gui.appdata.xml index abf0a4bf45..c019ab1ce3 100644 --- a/tools/lammps-gui/lammps-gui.appdata.xml +++ b/tools/lammps-gui/lammps-gui.appdata.xml @@ -61,6 +61,8 @@ Highlight warnings and error messages in Output window Make Tutorial wizards more compact Include download and compilation of WHAM software from Alan Grossfield + Add dialog to run WHAM directly from LAMMPS-GUI + Use mutex to avoid corruption of thermo data diff --git a/tools/lammps-gui/lammpsgui.cpp b/tools/lammps-gui/lammpsgui.cpp index 8749e24821..59e8017fc0 100644 --- a/tools/lammps-gui/lammpsgui.cpp +++ b/tools/lammps-gui/lammpsgui.cpp @@ -1015,6 +1015,7 @@ void LammpsGui::logupdate() void *ptr = lammps.last_thermo("setup", 0); if (ptr && *(int *)ptr) return; + lammps.last_thermo("lock", 0); ptr = lammps.last_thermo("num", 0); if (ptr) { int ncols = *(int *)ptr; @@ -1066,6 +1067,7 @@ void LammpsGui::logupdate() chartwindow->add_data(step, data, i); } } + lammps.last_thermo("unlock", 0); } // update list of available image file names From e20f3ec8745e81132ad71a07a221179ef8feee9d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 12 Jan 2025 23:50:33 -0500 Subject: [PATCH 31/37] macOS does not like forward declarations for standard C++ classes --- src/thermo.cpp | 1 - src/thermo.h | 5 +---- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/thermo.cpp b/src/thermo.cpp index d7a00fd14c..35b5016118 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -48,7 +48,6 @@ #include #include -#include #include using namespace LAMMPS_NS; diff --git a/src/thermo.h b/src/thermo.h index b5934f48b1..f5466e70cd 100644 --- a/src/thermo.h +++ b/src/thermo.h @@ -16,10 +16,7 @@ #include "pointers.h" #include - -namespace std { -class mutex; -} +#include namespace LAMMPS_NS { From a22c58cd47b2e5d80196e8a5e222e5a3bc7c9821 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Mon, 13 Jan 2025 00:17:05 -0700 Subject: [PATCH 32/37] add missing update to invoked_bonds in ComputeReaxFFAtomKokkos --- src/KOKKOS/compute_reaxff_atom_kokkos.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp index 0683e63752..ba6807641d 100644 --- a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp +++ b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp @@ -63,6 +63,8 @@ void ComputeReaxFFAtomKokkos::init() template void ComputeReaxFFAtomKokkos::compute_bonds() { + invoked_bonds = update->ntimestep; + if (atom->nmax > nmax) { memory->destroy(array_atom); nmax = atom->nmax; From 62aa803d53e6e30a93a9605a00a73801fe58db0a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 Jan 2025 12:03:31 -0500 Subject: [PATCH 33/37] order list of removed commands and packages in reverse order of their removal date also add an (automatic) table of contents simplify finding something --- doc/src/Commands_removed.rst | 182 +++++++++++++++++++---------------- 1 file changed, 100 insertions(+), 82 deletions(-) diff --git a/doc/src/Commands_removed.rst b/doc/src/Commands_removed.rst index 777be3bf16..be775f4c19 100644 --- a/doc/src/Commands_removed.rst +++ b/doc/src/Commands_removed.rst @@ -1,6 +1,10 @@ Removed commands and packages ============================= +.. contents:: + +------ + This page lists LAMMPS commands and packages that have been removed from the distribution and provides suggestions for alternatives or replacements. LAMMPS has special dummy styles implemented, that will @@ -8,47 +12,60 @@ stop LAMMPS and print a suitable error message in most cases, when a style/command is used that has been removed or will replace the command with the direct alternative (if available) and print a warning. -restart2data tool ------------------ +LAMMPS shell +------------ -.. versionchanged:: 23Nov2013 +.. versionchanged:: 29Aug2024 -The functionality of the restart2data tool has been folded into the -LAMMPS executable directly instead of having a separate tool. A -combination of the commands :doc:`read_restart ` and -:doc:`write_data ` can be used to the same effect. For -added convenience this conversion can also be triggered by -:doc:`command-line flags ` +The LAMMPS shell has been removed from the LAMMPS distribution. Users +are encouraged to use the :ref:`LAMMPS-GUI ` tool instead. -Fix ave/spatial and fix ave/spatial/sphere ------------------------------------------- +i-PI tool +--------- -.. deprecated:: 11Dec2015 +.. versionchanged:: 27Jun2024 -The fixes ave/spatial and ave/spatial/sphere have been removed from LAMMPS -since they were superseded by the more general and extensible "chunk -infrastructure". Here the system is partitioned in one of many possible -ways through the :doc:`compute chunk/atom ` command -and then averaging is done using :doc:`fix ave/chunk `. -Please refer to the :doc:`chunk HOWTO ` section for an overview. +The i-PI tool has been removed from the LAMMPS distribution. Instead, +instructions to install i-PI from PyPI via pip are provided. -Box command ------------ +USER-REAXC package +------------------ -.. deprecated:: 22Dec2022 +.. deprecated:: 7Feb2024 -The *box* command has been removed and the LAMMPS code changed so it won't -be needed. If present, LAMMPS will ignore the command and print a warning. +The USER-REAXC package has been renamed to :ref:`REAXFF `. +In the process also the pair style and related fixes were renamed to use +the "reaxff" string instead of "reax/c". For a while LAMMPS was maintaining +backward compatibility by providing aliases for the styles. These have +been removed, so using "reaxff" is now *required*. -Reset_ids, reset_atom_ids, reset_mol_ids commands -------------------------------------------------- +MPIIO package +------------- -.. deprecated:: 22Dec2022 +.. deprecated:: 21Nov2023 -The *reset_ids*, *reset_atom_ids*, and *reset_mol_ids* commands have -been folded into the :doc:`reset_atoms ` command. If -present, LAMMPS will replace the commands accordingly and print a -warning. +The MPIIO package has been removed from LAMMPS since it was unmaintained +for many years and thus not updated to incorporate required changes that +had been applied to the corresponding non-MPIIO commands. As a +consequence the MPIIO commands had become unreliable and sometimes +crashing LAMMPS or corrupting data. Similar functionality is available +through the :ref:`ADIOS package ` and the :ref:`NETCDF +package `. Also, the :doc:`dump_modify nfile or dump_modify +fileper ` keywords may be used for an efficient way of +writing out dump files when running on large numbers of processors. +Similarly, the "nfile" and "fileper" keywords exist for restarts: +see :doc:`restart `, :doc:`read_restart `, +:doc:`write_restart `. + +MSCG package +------------ + +.. deprecated:: 21Nov2023 + +The MSCG package has been removed from LAMMPS since it was unmaintained +for many years and instead superseded by the `OpenMSCG software +`_ of the Voth group at the +University of Chicago, which can be used independent from LAMMPS. LATTE package ------------- @@ -64,18 +81,6 @@ packages, including LATTE. See the ``examples/QUANTUM`` dir and the with LATTE as a plugin library (similar to the way fix latte worked), as well as on a different set of MPI processors. -MEAM package ------------- - -The MEAM package in Fortran has been replaced by a C++ implementation. -The code in the :ref:`MEAM package ` is a translation of the -Fortran code of MEAM into C++, which removes several restrictions -(e.g. there can be multiple instances in hybrid pair styles) and allows -for some optimizations leading to better performance. The pair style -:doc:`meam ` has the exact same syntax. For a transition -period the C++ version of MEAM was called USER-MEAMC so it could -coexist with the Fortran version. - Minimize style fire/old ----------------------- @@ -97,38 +102,38 @@ The same functionality is available through :doc:`bond style mesocnt ` and :doc:`angle style mesocnt `. -MPIIO package -------------- +Box command +----------- -.. deprecated:: 21Nov2023 +.. deprecated:: 22Dec2022 -The MPIIO package has been removed from LAMMPS since it was unmaintained -for many years and thus not updated to incorporate required changes that -had been applied to the corresponding non-MPIIO commands. As a -consequence the MPIIO commands had become unreliable and sometimes -crashing LAMMPS or corrupting data. Similar functionality is available -through the :ref:`ADIOS package ` and the :ref:`NETCDF -package `. Also, the :doc:`dump_modify nfile or dump_modify -fileper ` keywords may be used for an efficient way of -writing out dump files when running on large numbers of processors. -Similarly, the "nfile" and "fileper" keywords exist for restarts: -see :doc:`restart `, :doc:`read_restart `, -:doc:`write_restart `. +The *box* command has been removed and the LAMMPS code changed so it won't +be needed. If present, LAMMPS will ignore the command and print a warning. +Reset_ids, reset_atom_ids, reset_mol_ids commands +------------------------------------------------- -MSCG package ------------- +.. deprecated:: 22Dec2022 -.. deprecated:: 21Nov2023 +The *reset_ids*, *reset_atom_ids*, and *reset_mol_ids* commands have +been folded into the :doc:`reset_atoms ` command. If +present, LAMMPS will replace the commands accordingly and print a +warning. -The MSCG package has been removed from LAMMPS since it was unmaintained -for many years and instead superseded by the `OpenMSCG software -`_ of the Voth group at the -University of Chicago, which can be used independent from LAMMPS. +MESSAGE package +--------------- + +.. deprecated:: 4May2022 + +The MESSAGE package has been removed since it was superseded by the +:ref:`MDI package `. MDI implements the same functionality +and in a more general way with direct support for more applications. REAX package ------------ +.. deprecated:: 4Jan2019 + The REAX package has been removed since it was superseded by the :ref:`REAXFF package `. The REAXFF package has been tested to yield equivalent results to the REAX package, offers better @@ -138,20 +143,25 @@ syntax compatible with the removed reax pair style, so input files will have to be adapted. The REAXFF package was originally called USER-REAXC. -USER-REAXC package ------------------- +MEAM package +------------ -.. deprecated:: 7Feb2024 +.. deprecated:: 4Jan2019 -The USER-REAXC package has been renamed to :ref:`REAXFF `. -In the process also the pair style and related fixes were renamed to use -the "reaxff" string instead of "reax/c". For a while LAMMPS was maintaining -backward compatibility by providing aliases for the styles. These have -been removed, so using "reaxff" is now *required*. +The MEAM package in Fortran has been replaced by a C++ implementation. +The code in the :ref:`MEAM package ` is a translation of the +Fortran code of MEAM into C++, which removes several restrictions +(e.g. there can be multiple instances in hybrid pair styles) and allows +for some optimizations leading to better performance. The pair style +:doc:`meam ` has the exact same syntax. For a transition +period the C++ version of MEAM was called USER-MEAMC so it could +coexist with the Fortran version. USER-CUDA package ----------------- +.. deprecated:: 31May2016 + The USER-CUDA package had been removed, since it had been unmaintained for a long time and had known bugs and problems. Significant parts of the design were transferred to the @@ -160,19 +170,27 @@ performance characteristics on NVIDIA GPUs. Both, the KOKKOS and the :ref:`GPU package ` are maintained and allow running LAMMPS with GPU acceleration. -i-PI tool ---------- +Fix ave/spatial and fix ave/spatial/sphere +------------------------------------------ -.. versionchanged:: 27Jun2024 +.. deprecated:: 11Dec2015 -The i-PI tool has been removed from the LAMMPS distribution. Instead, -instructions to install i-PI from PyPI via pip are provided. +The fixes ave/spatial and ave/spatial/sphere have been removed from LAMMPS +since they were superseded by the more general and extensible "chunk +infrastructure". Here the system is partitioned in one of many possible +ways through the :doc:`compute chunk/atom ` command +and then averaging is done using :doc:`fix ave/chunk `. +Please refer to the :doc:`chunk HOWTO ` section for an overview. -LAMMPS shell ------------- +restart2data tool +----------------- -.. versionchanged:: 29Aug2024 +.. deprecated:: 23Nov2013 -The LAMMPS shell has been removed from the LAMMPS distribution. Users -are encouraged to use the :ref:`LAMMPS-GUI ` tool instead. +The functionality of the restart2data tool has been folded into the +LAMMPS executable directly instead of having a separate tool. A +combination of the commands :doc:`read_restart ` and +:doc:`write_data ` can be used to the same effect. For +added convenience this conversion can also be triggered by +:doc:`command-line flags ` From be048fc6360123a7b86fa14d4ead5f2e9a118776 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 Jan 2025 16:01:18 -0500 Subject: [PATCH 34/37] cosmetic changes --- doc/src/Library_objects.rst | 4 ++-- python/lammps/core.py | 2 +- python/lammps/ipython/wrapper.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/Library_objects.rst b/doc/src/Library_objects.rst index 89056a5db0..53edfce9e6 100644 --- a/doc/src/Library_objects.rst +++ b/doc/src/Library_objects.rst @@ -1,5 +1,5 @@ -Compute, fixes, variables -========================= +Computes, fixes, variables +========================== This section documents accessing or modifying data stored by computes, fixes, or variables in LAMMPS using the following functions: diff --git a/python/lammps/core.py b/python/lammps/core.py index accf20b5d1..3eb74cb550 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -103,7 +103,7 @@ class command_wrapper(object): This method is where the Python 'magic' happens. If a method is not defined by the class command_wrapper, it assumes it is a LAMMPS command. It takes all the arguments, concatinates them to a single string, and executes it using - :py:meth:`lammps.command()`. + :py:meth:`lammps.command`. Starting with Python 3.6 it also supports keyword arguments. key=value is transformed into 'key value'. Note, since these have come last in the diff --git a/python/lammps/ipython/wrapper.py b/python/lammps/ipython/wrapper.py index 729c0d62bf..6f6669f55a 100644 --- a/python/lammps/ipython/wrapper.py +++ b/python/lammps/ipython/wrapper.py @@ -17,7 +17,7 @@ ################################################################################ class wrapper(object): - """lammps API IPython Wrapper + """ lammps API IPython Wrapper This is a wrapper class that provides additional methods on top of an existing :py:class:`lammps` instance. It provides additional methods From 663f812799a7d8c97622dcec999d4643d0700d25 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 Jan 2025 19:03:05 -0500 Subject: [PATCH 35/37] fix inconsistent backquotes --- src/library.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/library.cpp b/src/library.cpp index aa835d0906..c5341eec83 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -517,7 +517,7 @@ The function returns the expanded string in a new string buffer that must be freed with :cpp:func:`lammps_free` after use to avoid a memory leak. *See also* - :cpp:func:`lammps_eval' + :cpp:func:`lammps_eval` \endverbatim * @@ -2936,7 +2936,7 @@ for :doc:`equal style variables `, evaluates it and returns the resulting (scalar) value as a floating point number. *See also* - :cpp:func:`lammps_expand' + :cpp:func:`lammps_expand` \endverbatim From 6d7926a026cb3eef5393be7a30d03910a7451db2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 Jan 2025 19:13:34 -0500 Subject: [PATCH 36/37] add missing lammps_expand() function --- doc/src/Library_execute.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/doc/src/Library_execute.rst b/doc/src/Library_execute.rst index 44b601ba4c..8560fbb148 100644 --- a/doc/src/Library_execute.rst +++ b/doc/src/Library_execute.rst @@ -7,6 +7,7 @@ This section documents the following functions: - :cpp:func:`lammps_command` - :cpp:func:`lammps_commands_list` - :cpp:func:`lammps_commands_string` +- :cpp:func:`lammps_expand` -------------------- @@ -79,3 +80,8 @@ Below is a short example using some of these functions. .. doxygenfunction:: lammps_commands_string :project: progguide +----------------------- + +.. doxygenfunction:: lammps_expand + :project: progguide + From 7c8c8c9d01b17d57ea661d0ffdfad87bb0ca533d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 13 Jan 2025 22:07:44 -0500 Subject: [PATCH 37/37] update patch for WHAM code --- tools/lammps-gui/update-wham-2.0.11.patch | 494 +++++++++++++++++++++- 1 file changed, 492 insertions(+), 2 deletions(-) diff --git a/tools/lammps-gui/update-wham-2.0.11.patch b/tools/lammps-gui/update-wham-2.0.11.patch index 3b6e0e372a..b884c25336 100644 --- a/tools/lammps-gui/update-wham-2.0.11.patch +++ b/tools/lammps-gui/update-wham-2.0.11.patch @@ -55,8 +55,79 @@ index 0000000..b4f0fe6 + TYPE DOC + PERMISSIONS OWNER_READ GROUP_READ WORLD_READ +) +diff --git a/nr/locate.c b/nr/locate.c +index 9f92dc0..f3bf294 100644 +--- a/nr/locate.c ++++ b/nr/locate.c +@@ -11,7 +11,7 @@ void locate(double xx[], int n, double x, int *j) + ascnd=(xx[n] > xx[0]); // I think this makes it zero based + while (ju-jl > 1) { + jm=(ju+jl) >> 1; +- if (x > xx[jm] == ascnd) ++ if ((x > xx[jm]) == ascnd) + jl=jm; + else + ju=jm; +diff --git a/wham-2d/histogram.c b/wham-2d/histogram.c +index 1bd1329..b5d1c01 100644 +--- a/wham-2d/histogram.c ++++ b/wham-2d/histogram.c +@@ -78,14 +78,14 @@ return hp; + } + + /* Get a value from a histogram structure +- * Given i and j, the indices into the global histogram, return ++ * Given i and j, the indices into the global histogram, return + * the correct histogram value. Since we only store range of histogram +- * indices containing the nonzero values, we have to check the index value ++ * indices containing the nonzero values, we have to check the index value + * against the range in the structure. + */ + double get_histval(struct histogram *hist, int i, int j) + { +-if ( (i < hist->first_x) || (i > hist->last_x) ++if ( (i < hist->first_x) || (i > hist->last_x) + ||(j < hist->first_y) || (j > hist->last_y)) + { + return 0.0; +@@ -239,7 +239,7 @@ return error; + + + // Calculate the free energy, setting the minimum value to 0 +-void calc_free(double **free, double **prob, double kT, ++void calc_free(double **free, double **prob, double kT, + int use_mask, int **mask) + { + int i,j; +@@ -257,7 +257,14 @@ for (i=0; i 0.0) ++ { ++ free[i][j] = -kT * log(prob[i][j]); ++ } ++ else ++ { ++ free[i][j] = 0.0; ++ } + if (free[i][j] < min) + { + min = free[i][j]; +@@ -321,8 +328,6 @@ return 0.5*(springx*(dx*dx) + springy*(dy*dy)); + + void calc_coor(int i, int j, double *coor) + { +-coor[0] = HIST_MINx + BIN_WIDTHx*((double)i+0.5); +-coor[1] = HIST_MINy + BIN_WIDTHy*((double)j+0.5); ++coor[0] = HIST_MINx + BIN_WIDTHx*((double)i+0.5); ++coor[1] = HIST_MINy + BIN_WIDTHy*((double)j+0.5); + } +- +- diff --git a/wham-2d/wham-2d.c b/wham-2d/wham-2d.c -index fb6e059..2c5594f 100644 +index fb6e059..a6b8483 100644 --- a/wham-2d/wham-2d.c +++ b/wham-2d/wham-2d.c @@ -25,7 +25,7 @@ @@ -77,6 +148,15 @@ index fb6e059..2c5594f 100644 int main(int argc, char *argv[]) { +@@ -57,7 +57,7 @@ double sum; + int iteration; + int max_iteration = 100000; + int numpad; +-int **mask; ++int **mask = NULL; + int use_mask; + + cpu1 = ((double) clock())/CLOCKS_PER_SEC; @@ -76,6 +76,61 @@ for (i=0; i 0.0) ++ { ++ free[i] = -kT * log(prob[i]); ++ } ++ else ++ { ++ free[i] = 0.0; ++ } ++ ++ if (free[i] < min) + { + min = free[i]; + bin_min = i; +@@ -252,5 +260,5 @@ return 0.5*dx*dx*spring; + + double calc_coor(int i) + { +-return HIST_MIN + BIN_WIDTH*((double)i+0.5); ++return HIST_MIN + BIN_WIDTH*((double)i+0.5); + } diff --git a/wham/wham.c b/wham/wham.c -index 487871b..1496eed 100644 +index 487871b..edb8125 100644 --- a/wham/wham.c +++ b/wham/wham.c +@@ -1,4 +1,4 @@ +-/* ++/* + * WHAM -- Weighted Histogram Analysis Method + * + * Reference 1: Computer Physics Communications, 91(1995) 275-282, Benoit Roux @@ -21,7 +21,7 @@ #include "wham.h" @@ -184,6 +458,15 @@ index 487871b..1496eed 100644 int main(int argc, char *argv[]) { +@@ -41,7 +42,7 @@ int first; + int bin_min; + int have_energy; + char *freefile; +-FILE *METAFILE, *FREEFILE; ++FILE *METAFILE, *FREEFILE; + struct hist_group *hist_group; + struct histogram *hp; + double coor; @@ -82,6 +83,61 @@ for (i=0; inum_windows); + // Figure out if we have trajectories at different temperatures. + // Missing temperatures are set to -1 in read_metadata, and + // since we require that either all trajectories specify a temperature +-// or all trajectories are assumed to be at the wham temperature, we only ++// or all trajectories are assumed to be at the wham temperature, we only + // have to check one of them + if (hist_group->kT[0] > 0) + { +@@ -257,7 +313,7 @@ if (hist_group->kT[0] > 0) + else + { + have_energy = 0; +- for (i=0; i< hist_group->num_windows; i++) ++ for (i=0; i< hist_group->num_windows; i++) + { + hist_group->kT[i] = kT; + } +@@ -269,7 +325,7 @@ if (!final_f) + { + printf("couldn't allocate space for final_f: %s\n", strerror(errno)); + exit(errno); +- } ++ } + + free(HISTOGRAM); + +@@ -305,7 +361,8 @@ while (! is_converged(hist_group) || first ) + for (i=0; i< NUM_BINS; i++) + { + coor = calc_coor(i); +- printf("%f\t%f\t%f\n", coor, free_ene[i], prob[i]); ++ if (prob[i] != 0.0) ++ printf("%f\t%f\t%f\n", coor, free_ene[i], prob[i]); + } + printf("\n"); + +@@ -319,7 +376,7 @@ while (! is_converged(hist_group) || first ) + } + } + // Cheesy bailout if we're going on too long +- if (iteration >= max_iteration) ++ if (iteration >= max_iteration) + { + printf("Too many iterations: %d\n", iteration); + break; +@@ -383,11 +440,11 @@ for (i=0; i< num_mc_runs; i++) + //printf("Faking %d: %d %d\n", i,j,hp->num_points); + num_used = hp->last - hp->first + 1; + mk_new_hist(hp->cum, hp->data, num_used, hp->num_mc_samples, &idum); +- ++ + hist_group->F_old[j] = 0.0; + hist_group->F[j] = 0.0; + } +- ++ + // perform WHAM iterations on the fake data sets + iteration = 0; + first = 1; +@@ -403,7 +460,7 @@ for (i=0; i< num_mc_runs; i++) + printf("Too many iterations: %d\n", iteration); + break; + } +- } ++ } + printf("#MC trial %d: %d iterations\n", i, iteration); + printf("#PMF values\n"); + // accumulate the average and stdev of the resulting probabilities +@@ -419,18 +476,19 @@ for (i=0; i< num_mc_runs; i++) + for (j=0; j < NUM_BINS; j++) + { + pdf = -kT*log(prob[j]); +- ++ + ave_p[j] += prob[j]; + ave_pdf[j] += pdf; + ave_p2[j] += prob[j] * prob[j]; + ave_pdf2[j] += pdf*pdf; + } +- for (j=0; jnum_windows;j++) ++ for (j=0; jnum_windows;j++) + { + ave_F[j] += hist_group->F[j] - hist_group->F[0]; +- ave_F2[j] += hist_group->F[j]*hist_group->F[j] ; ++ ave_F2[j] += hist_group->F[j]*hist_group->F[j] ; + } +- } ++ } ++ if (num_mc_runs == 0) num_mc_runs = 1; + for (i=0; i < NUM_BINS; i++) + { + ave_p[i] /= (double)num_mc_runs; +@@ -457,12 +515,12 @@ if (!FREEFILE) + for (i=0; i< NUM_BINS; i++) + { + coor = calc_coor(i); +- printf("%f\t%f\t%f\t%f\t%f\n", coor, free_ene[i], ave_pdf2[i], +- prob[i], ave_p2[i]); ++ if (prob[i] != 0.0) ++ printf("%f\t%f\t%f\t%f\t%f\n", coor, free_ene[i], ave_pdf2[i], prob[i], ave_p2[i]); + } + for (i=0; inum_windows; i++) + { +- fprintf(FREEFILE,"%d\t%f\t%f\n", i, final_f[i],ave_F2[i]); ++ printf("%d\t%f\t%f\n", i, final_f[i],ave_F2[i]); + } + + exit(errno); +@@ -470,38 +528,37 @@ if (!FREEFILE) + else + { + // write out header +- fprintf(FREEFILE, "#Coor\t\tFree\t+/-\t\tProb\t\t+/-\n"); ++ fprintf(FREEFILE, "#Coor\t\tFree\t\t+/-\t\tProb\t\t+/-\n"); + // write out the leading padded values + for (i=-numpad; i<0; i++) + { + coor = calc_coor(i); +- fprintf(FREEFILE,"%f\t%f\t%f\t%f\t%f\n", coor, free_ene[NUM_BINS+i], +- ave_pdf2[NUM_BINS+i], +- final_prob[NUM_BINS+i], +- ave_p2[NUM_BINS+i]); ++ if (prob[NUM_BINS+1] != 0.0) ++ fprintf(FREEFILE,"%f\t%f\t%f\t%f\t%f\n", coor, free_ene[NUM_BINS+i], ++ ave_pdf2[NUM_BINS+i], final_prob[NUM_BINS+i], ave_p2[NUM_BINS+i]); + } + // write out the center values + for (i=0; inum_windows; i++) + { +- fprintf(FREEFILE,"#%d\t%f\t%f\n", i, final_f[i],ave_F2[i]); ++ fprintf(FREEFILE,"#%d\t\t%f\t%f\n", i, final_f[i],ave_F2[i]); + } + } + +@@ -515,7 +572,7 @@ exit(0); + /* + * Perform a single WHAM iteration + */ +-void wham_iteration(struct hist_group* hist_group, double *prob, ++void wham_iteration(struct hist_group* hist_group, double *prob, + int have_energy) + { + int i,j; +@@ -535,9 +592,9 @@ for (i=0; ihists[j]),i); + bias = calc_bias(hist_group,j,coor); + bf = exp((hist_group->F_old[j] - bias) / hist_group->kT[j]); +- /* If we have energies, then the histograms contain boltzmann +- * factors. If not, they contain integer counts. Accordingly, +- * we either use a partition function to normalize, or simply the ++ /* If we have energies, then the histograms contain boltzmann ++ * factors. If not, they contain integer counts. Accordingly, ++ * we either use a partition function to normalize, or simply the + * number of points. + */ + if (have_energy) diff --git a/wham/wham.h b/wham/wham.h index aacc1e8..7d509f2 100644 --- a/wham/wham.h