From 99c9091025b79492308436b7a673d7d68d7db932 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 17 Aug 2011 21:56:18 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6713 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-MISC/fix_imd.cpp | 1476 ++++++++++++++++++++++++++++++++++ src/USER-MISC/fix_imd.h | 131 +++ src/USER-MISC/fix_smd.cpp | 403 ++++++++++ src/USER-MISC/fix_smd.h | 60 ++ src/USER-MISC/pair_cdeam.cpp | 642 +++++++++++++++ src/USER-MISC/pair_cdeam.h | 230 ++++++ 6 files changed, 2942 insertions(+) create mode 100644 src/USER-MISC/fix_imd.cpp create mode 100644 src/USER-MISC/fix_imd.h create mode 100644 src/USER-MISC/fix_smd.cpp create mode 100644 src/USER-MISC/fix_smd.h create mode 100644 src/USER-MISC/pair_cdeam.cpp create mode 100644 src/USER-MISC/pair_cdeam.h diff --git a/src/USER-MISC/fix_imd.cpp b/src/USER-MISC/fix_imd.cpp new file mode 100644 index 0000000000..3b34ac69b6 --- /dev/null +++ b/src/USER-MISC/fix_imd.cpp @@ -0,0 +1,1476 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + The FixIMD class contains code from VMD and NAMD which is copyrighted + by the Board of Trustees of the University of Illinois and is free to + use with LAMMPS according to point 2 of the UIUC license (10% clause): + +" Licensee may, at its own expense, create and freely distribute +complimentary works that interoperate with the Software, directing others to +the TCBG server to license and obtain the Software itself. Licensee may, at +its own expense, modify the Software to make derivative works. Except as +explicitly provided below, this License shall apply to any derivative work +as it does to the original Software distributed by Illinois. Any derivative +work should be clearly marked and renamed to notify users that it is a +modified version and not the original Software distributed by Illinois. +Licensee agrees to reproduce the copyright notice and other proprietary +markings on any derivative work and to include in the documentation of such +work the acknowledgement: + + "This software includes code developed by the Theoretical and Computational + Biophysics Group in the Beckman Institute for Advanced Science and + Technology at the University of Illinois at Urbana-Champaign." + +Licensee may redistribute without restriction works with up to 1/2 of their +non-comment source code derived from at most 1/10 of the non-comment source +code developed by Illinois and contained in the Software, provided that the +above directions for notice and acknowledgement are observed. Any other +distribution of the Software or any derivative work requires a separate +license with Illinois. Licensee may contact Illinois (vmd@ks.uiuc.edu) to +negotiate an appropriate license for such distribution." +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (TempleU) + IMD API, hash, and socket code written by: John E. Stone, + Justin Gullingsrud, and James Phillips, (TCBG, Beckman Institute, UIUC) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include + +#if defined(_MSC_VER) || defined(__MINGW32_VERSION) +#include +#else +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +#include + +#include "fix_imd.h" +#include "atom.h" +#include "comm.h" +#include "update.h" +#include "respa.h" +#include "domain.h" +#include "error.h" +#include "group.h" +#include "memory.h" + +/* re-usable integer hash table code with static linkage. */ + +/** hash table top level data structure */ +typedef struct inthash_t { + struct inthash_node_t **bucket; /* array of hash nodes */ + int size; /* size of the array */ + int entries; /* number of entries in table */ + int downshift; /* shift cound, used in hash function */ + int mask; /* used to select bits for hashing */ +} inthash_t; + +/** hash table node data structure */ +typedef struct inthash_node_t { + int data; /* data in hash node */ + int key; /* key for hash lookup */ + struct inthash_node_t *next; /* next node in hash chain */ +} inthash_node_t; + +#define HASH_FAIL -1 +#define HASH_LIMIT 0.5 + +/* initialize new hash table */ +static void inthash_init(inthash_t *tptr, int buckets); +/* lookup entry in hash table */ +static int inthash_lookup(const inthash_t *tptr, int key); +/* generate list of keys for reverse lookups. */ +static int *inthash_keys(inthash_t *tptr); +/* insert an entry into hash table. */ +static int inthash_insert(inthash_t *tptr, int key, int data); +/* delete the hash table */ +static void inthash_destroy(inthash_t *tptr); +/* adapted sort for in-place sorting of map indices. */ +static void id_sort(int *idmap, int left, int right); + +/************************************************************************ + * integer hash code: + ************************************************************************/ + +/* inthash() - Hash function returns a hash number for a given key. + * tptr: Pointer to a hash table, key: The key to create a hash number for */ +static int inthash(const inthash_t *tptr, int key) { + int hashvalue; + + hashvalue = (((key*1103515249)>>tptr->downshift) & tptr->mask); + if (hashvalue < 0) { + hashvalue = 0; + } + + return hashvalue; +} + +/* + * rebuild_table_int() - Create new hash table when old one fills up. + * + * tptr: Pointer to a hash table + */ +static void rebuild_table_int(inthash_t *tptr) { + inthash_node_t **old_bucket, *old_hash, *tmp; + int old_size, h, i; + + old_bucket=tptr->bucket; + old_size=tptr->size; + + /* create a new table and rehash old buckets */ + inthash_init(tptr, old_size<<1); + for (i=0; inext; + h=inthash(tptr, tmp->key); + tmp->next=tptr->bucket[h]; + tptr->bucket[h]=tmp; + tptr->entries++; + } /* while */ + } /* for */ + + /* free memory used by old table */ + free(old_bucket); + + return; +} + +/* + * inthash_init() - Initialize a new hash table. + * + * tptr: Pointer to the hash table to initialize + * buckets: The number of initial buckets to create + */ +void inthash_init(inthash_t *tptr, int buckets) { + + /* make sure we allocate something */ + if (buckets==0) + buckets=16; + + /* initialize the table */ + tptr->entries=0; + tptr->size=2; + tptr->mask=1; + tptr->downshift=29; + + /* ensure buckets is a power of 2 */ + while (tptr->sizesize<<=1; + tptr->mask=(tptr->mask<<1)+1; + tptr->downshift--; + } /* while */ + + /* allocate memory for table */ + tptr->bucket=(inthash_node_t **) calloc(tptr->size, sizeof(inthash_node_t *)); + + return; +} + +/* + * inthash_lookup() - Lookup an entry in the hash table and return a pointer to + * it or HASH_FAIL if it wasn't found. + * + * tptr: Pointer to the hash table + * key: The key to lookup + */ +int inthash_lookup(const inthash_t *tptr, int key) { + int h; + inthash_node_t *node; + + + /* find the entry in the hash table */ + h=inthash(tptr, key); + for (node=tptr->bucket[h]; node!=NULL; node=node->next) { + if (node->key == key) + break; + } + + /* return the entry if it exists, or HASH_FAIL */ + return(node ? node->data : HASH_FAIL); +} + + +/* + * inthash_keys() - Return a list of keys. + * NOTE: the returned list must be freed with free(3). + */ +int *inthash_keys(inthash_t *tptr) { + + int *keys; + inthash_node_t *node; + + keys = (int *)calloc(tptr->entries, sizeof(int)); + + for (int i=0; i < tptr->size; ++i) { + for (node=tptr->bucket[i]; node != NULL; node=node->next) { + keys[node->data] = node->key; + } + } + + return keys; +} + +/* + * inthash_insert() - Insert an entry into the hash table. If the entry already + * exists return a pointer to it, otherwise return HASH_FAIL. + * + * tptr: A pointer to the hash table + * key: The key to insert into the hash table + * data: A pointer to the data to insert into the hash table + */ +int inthash_insert(inthash_t *tptr, int key, int data) { + int tmp; + inthash_node_t *node; + int h; + + /* check to see if the entry exists */ + if ((tmp=inthash_lookup(tptr, key)) != HASH_FAIL) + return(tmp); + + /* expand the table if needed */ + while (tptr->entries>=HASH_LIMIT*tptr->size) + rebuild_table_int(tptr); + + /* insert the new entry */ + h=inthash(tptr, key); + node=(struct inthash_node_t *) malloc(sizeof(inthash_node_t)); + node->data=data; + node->key=key; + node->next=tptr->bucket[h]; + tptr->bucket[h]=node; + tptr->entries++; + + return HASH_FAIL; +} + +/* + * inthash_destroy() - Delete the entire table, and all remaining entries. + * + */ +void inthash_destroy(inthash_t *tptr) { + inthash_node_t *node, *last; + int i; + + for (i=0; isize; i++) { + node = tptr->bucket[i]; + while (node != NULL) { + last = node; + node = node->next; + free(last); + } + } + + /* free the entire array of buckets */ + if (tptr->bucket != NULL) { + free(tptr->bucket); + memset(tptr, 0, sizeof(inthash_t)); + } +} + +/************************************************************************ + * integer list sort code: + ************************************************************************/ + +/* sort for integer map. initial call id_sort(idmap, 0, natoms - 1); */ +static void id_sort(int *idmap, int left, int right) +{ + int pivot, l_hold, r_hold; + + l_hold = left; + r_hold = right; + pivot = idmap[left]; + + while (left < right) { + while ((idmap[right] >= pivot) && (left < right)) + right--; + if (left != right) { + idmap[left] = idmap[right]; + left++; + } + while ((idmap[left] <= pivot) && (left < right)) + left++; + if (left != right) { + idmap[right] = idmap[left]; + right--; + } + } + idmap[left] = pivot; + pivot = left; + left = l_hold; + right = r_hold; + + if (left < pivot) + id_sort(idmap, left, pivot-1); + if (right > pivot) + id_sort(idmap, pivot+1, right); +} + +/********** API definitions of the VMD/NAMD code ************************ + * This code was taken and adapted from VMD-1.8.7/NAMD-2.7 in Sep 2009. * + * If there are any bugs or problems, please contact akohlmey@gmail.com * + ************************************************************************/ + +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2009 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ + +/* part 1: Interactive MD (IMD) API */ + +#include + +#if ( INT_MAX == 2147483647 ) +typedef int int32; +#else +typedef short int32; +#endif + +typedef struct { + int32 type; + int32 length; +} IMDheader; + +#define IMDHEADERSIZE 8 +#define IMDVERSION 2 + +typedef enum IMDType_t { + IMD_DISCONNECT, /**< close IMD connection, leaving sim running */ + IMD_ENERGIES, /**< energy data block */ + IMD_FCOORDS, /**< atom coordinates */ + IMD_GO, /**< start the simulation */ + IMD_HANDSHAKE, /**< endianism and version check message */ + IMD_KILL, /**< kill the simulation job, shutdown IMD */ + 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 */ +} IMDType; /**< IMD command message type enumerations */ + +typedef struct { + int32 tstep; /**< integer timestep index */ + float T; /**< Temperature in degrees Kelvin */ + float Etot; /**< Total energy, in Kcal/mol */ + float Epot; /**< Potential energy, in Kcal/mol */ + float Evdw; /**< Van der Waals energy, in Kcal/mol */ + float Eelec; /**< Electrostatic energy, in Kcal/mol */ + float Ebond; /**< Bond energy, Kcal/mol */ + float Eangle; /**< Angle energy, Kcal/mol */ + float Edihe; /**< Dihedral energy, Kcal/mol */ + float Eimpr; /**< Improper energy, Kcal/mol */ +} IMDEnergies; /**< IMD simulation energy report structure */ + +/** Send control messages - these consist of a header with no subsequent data */ +static int imd_handshake(void *); /**< check endianness, version compat */ +/** Receive header and data */ +static IMDType imd_recv_header(void *, int32 *); +/** Receive MDComm-style forces, units are Kcal/mol/angstrom */ +static int imd_recv_mdcomm(void *, int32, int32 *, float *); +/** Receive energies */ +static int imd_recv_energies(void *, IMDEnergies *); +/** Receive atom coordinates. */ +static int imd_recv_fcoords(void *, int32, float *); +/** Prepare IMD data packet header */ +static void imd_fill_header(IMDheader *header, IMDType type, int32 length); +/** Write data to socket */ +static int32 imd_writen(void *s, const char *ptr, int32 n); + +/* part 2: abstracts platform-dependent routines/APIs for using sockets */ + +typedef struct { + struct sockaddr_in addr; /* address of socket provided by bind() */ + int addrlen; /* size of the addr struct */ + int sd; /* socket file descriptor */ +} imdsocket; + +static int imdsock_init(void); +static void *imdsock_create(void); +static int imdsock_bind(void *, int); +static int imdsock_listen(void *); +static void *imdsock_accept(void *); /* return new socket */ +static int imdsock_write(void *, const void *, int); +static int imdsock_read(void *, void *, int); +static int imdsock_selread(void *, int); +static int imdsock_selwrite(void *, int); +static void imdsock_shutdown(void *); +static void imdsock_destroy(void *); + +/*************************************************************** + * End of API definitions of the VMD/NAMD code. * + * The implementation follows at the end of the file. * + ***************************************************************/ + +using namespace LAMMPS_NS; + +/* struct for packed data communication of coordinates and forces. */ +struct commdata { + int tag; + float x,y,z; +}; + +/*************************************************************** + * create class and parse arguments in LAMMPS script. Syntax: + * fix ID group-ID imd [unwrap (on|off)] [fscale ] + ***************************************************************/ +FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 4) + error->all("Illegal fix imd command"); + + imd_port = atoi(arg[3]); + if (imd_port < 1024) + error->all("Illegal fix imd parameter: port < 1024"); + + /* default values for optional flags */ + unwrap_flag = 0; + nowait_flag = 0; + connect_msg = 1; + imd_fscale = 1.0; + imd_trate = 1; + + /* parse optional arguments */ + int argsdone = 4; + while (argsdone+1 < narg) { + if (0 == strcmp(arg[argsdone], "unwrap")) { + if (0 == strcmp(arg[argsdone+1], "on")) { + unwrap_flag = 1; + } else { + unwrap_flag = 0; + } + } else if (0 == strcmp(arg[argsdone], "nowait")) { + if (0 == strcmp(arg[argsdone+1], "on")) { + nowait_flag = 1; + } else { + nowait_flag = 0; + } + } else if (0 == strcmp(arg[argsdone], "fscale")) { + imd_fscale = atof(arg[argsdone+1]); + } else if (0 == strcmp(arg[argsdone], "trate")) { + imd_trate = atoi(arg[argsdone+1]); + } else { + error->all("Unknown fix imd parameter"); + } + ++argsdone; ++argsdone; + } + + /* sanity check on parameters */ + if (imd_trate < 1) + error->all("Illegal fix imd parameter. trate < 1."); + + bigint n = group->count(igroup); + if (n > MAXSMALLINT) error->all("Too many atoms for fix imd"); + num_coords = static_cast (n); + + MPI_Comm_rank(world,&me); + + /* initialize various imd state variables. */ + clientsock = NULL; + localsock = NULL; + nlevels_respa = 0; + imd_inactive = 0; + imd_terminate = 0; + imd_forces = 0; + force_buf = NULL; + maxbuf = 0; + msgdata = NULL; + msglen = 0; + comm_buf = NULL; + idmap = NULL; + rev_idmap = NULL; + + if (me == 0) { + /* set up incoming socket on MPI rank 0. */ + imdsock_init(); + localsock = imdsock_create(); + clientsock = NULL; + if (imdsock_bind(localsock,imd_port)) { + perror("bind to socket failed"); + imdsock_destroy(localsock); + imd_terminate = 1; + } else { + imdsock_listen(localsock); + } + } + MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); + if (imd_terminate) + error->all("LAMMPS Terminated on error in IMD."); + + /* storage required to communicate a single coordinate or force. */ + size_one = sizeof(struct commdata); + +#if defined(LAMMPS_ASYNC_IMD) + /* set up for i/o worker thread on MPI rank 0.*/ + if (me == 0) { + if (screen) + fputs("Using fix imd with asynchronous I/O.\n",screen); + if (logfile) + fputs("Using fix imd with asynchronous I/O.\n",logfile); + + /* set up mutex and condition variable for i/o thread */ + /* hold mutex before creating i/o thread to keep it waiting. */ + pthread_mutex_init(&read_mutex, NULL); + pthread_mutex_init(&write_mutex, NULL); + pthread_cond_init(&write_cond, NULL); + + pthread_mutex_lock(&write_mutex); + buf_has_data=0; + pthread_mutex_unlock(&write_mutex); + + /* set up and launch i/o thread */ + pthread_attr_init(&iot_attr); + pthread_attr_setdetachstate(&iot_attr, PTHREAD_CREATE_JOINABLE); + pthread_create(&iothread, &iot_attr, &fix_imd_ioworker, this); + } +#endif +} + +/********************************* + * Clean up on deleting the fix. * + *********************************/ +FixIMD::~FixIMD() +{ + +#if defined(LAMMPS_ASYNC_IMD) + if (me == 0) { + pthread_mutex_lock(&write_mutex); + buf_has_data=-1; + pthread_cond_signal(&write_cond); + pthread_mutex_unlock(&write_mutex); + pthread_join(iothread, NULL); + + /* cleanup */ + pthread_attr_destroy(&iot_attr); + pthread_mutex_destroy(&write_mutex); + pthread_cond_destroy(&write_cond); + } +#endif + + inthash_t *hashtable = (inthash_t *)idmap; + memory->sfree(comm_buf); + memory->sfree(force_buf); + inthash_destroy(hashtable); + delete hashtable; + free(rev_idmap); + // close sockets + imdsock_shutdown(clientsock); + imdsock_destroy(clientsock); + imdsock_shutdown(localsock); + imdsock_destroy(localsock); + clientsock=NULL; + localsock=NULL; + return; +} + +/* ---------------------------------------------------------------------- */ +int FixIMD::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= POST_FORCE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ +void FixIMD::init() +{ + if (strstr(update->integrate_style,"respa")) + nlevels_respa = ((Respa *) update->integrate)->nlevels; + + return; +} + +/* ---------------------------------------------------------------------- */ + +/* (re-)connect to an IMD client (e.g. VMD). return 1 if + new connection was made, 0 if not. */ +int FixIMD::reconnect() +{ + /* set up IMD communication, but only if needed. */ + imd_inactive = 0; + imd_terminate = 0; + + if (me == 0) { + if (clientsock) return 1; + if (screen && connect_msg) + if (nowait_flag) + fprintf(screen,"Listening for IMD connection on port %d. Transfer rate %d.\n",imd_port, imd_trate); + else + fprintf(screen,"Waiting for IMD connection on port %d. Transfer rate %d.\n",imd_port, imd_trate); + + connect_msg = 0; + clientsock = NULL; + if (nowait_flag) { + int retval = imdsock_selread(localsock,0); + if (retval > 0) { + clientsock = imdsock_accept(localsock); + } else { + imd_inactive = 1; + return 0; + } + } else { + int retval=0; + do { + retval = imdsock_selread(localsock, 60); + } while (retval <= 0); + clientsock = imdsock_accept(localsock); + } + + if (!imd_inactive && !clientsock) { + if (screen) + fprintf(screen, "IMD socket accept error. Dropping connection.\n"); + imd_terminate = 1; + return 0; + } else { + /* check endianness and IMD protocol version. */ + if (imd_handshake(clientsock)) { + if (screen) + fprintf(screen, "IMD handshake error. Dropping connection.\n"); + imdsock_destroy(clientsock); + imd_terminate = 1; + return 0; + } else { + int32 length; + if (imdsock_selread(clientsock, 1) != 1 || + imd_recv_header(clientsock, &length) != IMD_GO) { + if (screen) + fprintf(screen, "Incompatible IMD client version? Dropping connection.\n"); + imdsock_destroy(clientsock); + imd_terminate = 1; + return 0; + } else { + return 1; + } + } + } + } + return 0; +} + +/* ---------------------------------------------------------------------- */ +/* wait for IMD client (e.g. VMD) to respond, initialize communication + * buffers and collect tag/id maps. */ +void FixIMD::setup(int) +{ + /* 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; + int *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); + maxbuf = nmax*size_one; + comm_buf = (void *) memory->smalloc(maxbuf,"imd:comm_buf"); + + 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("LAMMPS terminated on error in setting up IMD connection."); + + /* initialize and build hashtable. */ + inthash_t *hashtable=new inthash_t; + inthash_init(hashtable, num_coords); + idmap = (void *)hashtable; + + MPI_Status status; + MPI_Request request; + int tmp, ndata; + struct commdata *buf = static_cast(comm_buf); + + if (me == 0) { + int *taglist = new int[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(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; + + 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) { + inthash_insert(hashtable, taglist[i], i); + } + delete[] taglist; + + /* generate reverse index-to-tag map for communicating + * IMD forces back to the proper atoms */ + rev_idmap=inthash_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, &status); + MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world); + } + + return; +} + +/* worker threads for asynchronous i/o */ +#if defined(LAMMPS_ASYNC_IMD) +/* c bindings wrapper */ +void *fix_imd_ioworker(void *t) +{ + FixIMD *imd=(FixIMD *)t; + imd->ioworker(); + return NULL; +} + +/* the real i/o worker thread */ +void FixIMD::ioworker() +{ + while (1) { + pthread_mutex_lock(&write_mutex); + if (buf_has_data < 0) { + /* master told us to go away */ + fprintf(screen,"Asynchronous I/O thread is exiting.\n"); + buf_has_data=0; + pthread_mutex_unlock(&write_mutex); + pthread_exit(NULL); + } 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); + } + delete[] msgdata; + buf_has_data=0; + pthread_mutex_unlock(&write_mutex); + } else { + /* nothing to write out yet. wait on condition. */ + pthread_cond_wait(&write_cond, &write_mutex); + pthread_mutex_unlock(&write_mutex); + } + } +} +#endif + +/* ---------------------------------------------------------------------- */ +/* Main IMD protocol handler: + * Send coodinates, energies, and add IMD forces to atoms. */ +void FixIMD::post_force(int vflag) +{ + /* 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("LAMMPS terminated on error in setting up IMD connection."); + if (imd_inactive) + return; /* IMD client has detached and not yet come back. do nothing. */ + } + + int *tag = atom->tag; + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + int *mask = atom->mask; + struct commdata *buf; + + 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_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 = NULL; + imdsock_destroy(clientsock); + clientsock = NULL; + 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 = NULL; + break; + + 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; + } else { + if (screen) + fprintf(screen, "Pausing run on IMD client request.\n"); + imd_paused = 1; + } + 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_ENERGIES: { + IMDEnergies dummy_energies; + imd_recv_energies(clientsock, &dummy_energies); + break; + } + + case IMD_FCOORDS: { + float *dummy_coords = new float[3*length]; + imd_recv_fcoords(clientsock, length, dummy_coords); + delete[] dummy_coords; + break; + } + + case IMD_MDCOMM: { + int32 *imd_tags = new int32[length]; + float *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(length*size_one, + "imd:force_buf"); + } + imd_forces = length; + buf = static_cast(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; + } + + 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("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 (force_buf != NULL) + memory->sfree(force_buf); + force_buf = memory->smalloc(imd_forces*size_one, "imd:force_buf"); + } + } + MPI_Bcast(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 (update->ntimestep % imd_trate) { + if (imd_forces > 0) { + double **f = atom->f; + buf = static_cast(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; + } + } + } + } + } + return; + } + + /* check and potentially grow local communication buffers. */ + int i, k, nmax, nme=0; + for (i=0; i < nlocal; ++i) + if (mask[i] & groupbit) ++nme; + + MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); + if (nmax*size_one > maxbuf) { + memory->sfree(comm_buf); + maxbuf = nmax*size_one; + comm_buf = memory->smalloc(maxbuf,"imd:comm_buf"); + } + + MPI_Status status; + MPI_Request request; + int tmp, ndata; + buf = static_cast(comm_buf); + + if (me == 0) { + /* 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. */ + float *recvcoord = (float *) (msgdata+IMDHEADERSIZE); + + /* add local data */ + if (unwrap_flag) { + 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 (i=0; i> 10 & 1023) - 512; + int iz = (image[i] >> 20) - 512; + + if (domain->triclinic) { + recvcoord[j] = x[i][0] + ix * xprd + iy * xy + iz * xz; + recvcoord[j+1] = x[i][1] + iy * yprd + iz * yz; + recvcoord[j+2] = x[i][2] + iz * zprd; + } else { + recvcoord[j] = x[i][0] + ix * xprd; + recvcoord[j+1] = x[i][1] + iy * yprd; + recvcoord[j+2] = x[i][2] + iz * zprd; + } + } + } + } + } else { + for (i=0; inprocs; ++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; + + for (k=0; kxprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + double xy = domain->xy; + double xz = domain->xz; + double yz = domain->yz; + + for (i=0; i> 10 & 1023) - 512; + int iz = (image[i] >> 20) - 512; + + if (domain->triclinic) { + buf[nme].tag = tag[i]; + buf[nme].x = x[i][0] + ix * xprd + iy * xy + iz * xz; + buf[nme].y = x[i][1] + iy * yprd + iz * yz; + buf[nme].z = x[i][2] + iz * zprd; + } else { + buf[nme].tag = tag[i]; + buf[nme].x = x[i][0] + ix * xprd; + buf[nme].y = x[i][1] + iy * yprd; + buf[nme].z = x[i][2] + iz * zprd; + } + ++nme; + } + } + } else { + for (i=0; i(num_coords+maxbuf+imd_forces)*size_one; +} + +/* End of FixIMD class implementation. */ + +/***************************************************************************/ + +/* NOTE: the following code is the based on the example implementation + * of the IMD protocol API from VMD and NAMD. The UIUC license allows + * to re-use up to 10% of a project's code to be used in other software */ + +/*************************************************************************** + * DESCRIPTION: + * Socket interface, abstracts machine dependent APIs/routines. + ***************************************************************************/ + +int imdsock_init(void) { +#if defined(_MSC_VER) || defined(__MINGW32_VERSION) + int rc = 0; + static int initialized=0; + + if (!initialized) { + WSADATA wsdata; + rc = WSAStartup(MAKEWORD(1,1), &wsdata); + if (rc == 0) + initialized = 1; + } + + return rc; +#else + return 0; +#endif +} + + +void * imdsock_create(void) { + imdsocket * s; + + s = (imdsocket *) malloc(sizeof(imdsocket)); + if (s != NULL) + memset(s, 0, sizeof(imdsocket)); + + if ((s->sd = socket(PF_INET, SOCK_STREAM, 0)) == -1) { + printf("Failed to open socket."); + free(s); + return NULL; + } + + return (void *) s; +} + +int imdsock_bind(void * v, int port) { + imdsocket *s = (imdsocket *) v; + memset(&(s->addr), 0, sizeof(s->addr)); + s->addr.sin_family = PF_INET; + s->addr.sin_port = htons(port); + + return bind(s->sd, (struct sockaddr *) &s->addr, sizeof(s->addr)); +} + +int imdsock_listen(void * v) { + imdsocket *s = (imdsocket *) v; + return listen(s->sd, 5); +} + +void *imdsock_accept(void * v) { + int rc; + imdsocket *new_s = NULL, *s = (imdsocket *) v; +#if defined(ARCH_AIX5) || defined(ARCH_AIX5_64) || defined(ARCH_AIX6_64) + unsigned int len; +#define _SOCKLEN_TYPE unsigned int +#elif defined(SOCKLEN_T) + SOCKLEN_T len; +#define _SOCKLEN_TYPE SOCKLEN_T +#elif defined(_POSIX_SOURCE) + socklen_t len; +#define _SOCKLEN_TYPE socklen_t +#else +#define _SOCKLEN_TYPE int + int len; +#endif + + len = sizeof(s->addr); + rc = accept(s->sd, (struct sockaddr *) &s->addr, ( _SOCKLEN_TYPE * ) &len); + if (rc >= 0) { + new_s = (imdsocket *) malloc(sizeof(imdsocket)); + if (new_s != NULL) { + *new_s = *s; + new_s->sd = rc; + } + } + return (void *)new_s; +} + +int imdsock_write(void * v, const void *buf, int len) { + imdsocket *s = (imdsocket *) v; +#if defined(_MSC_VER) || defined(__MINGW32_VERSION) + return send(s->sd, (const char*) buf, len, 0); /* windows lacks the write() call */ +#else + return write(s->sd, buf, len); +#endif +} + +int imdsock_read(void * v, void *buf, int len) { + imdsocket *s = (imdsocket *) v; +#if defined(_MSC_VER) || defined(__MINGW32_VERSION) + return recv(s->sd, (char*) buf, len, 0); /* windows lacks the read() call */ +#else + return read(s->sd, buf, len); +#endif + +} + +void imdsock_shutdown(void *v) { + imdsocket * s = (imdsocket *) v; + if (s == NULL) + return; + +#if defined(_MSC_VER) || defined(__MINGW32_VERSION) + shutdown(s->sd, SD_SEND); +#else + shutdown(s->sd, 1); /* complete sends and send FIN */ +#endif +} + +void imdsock_destroy(void * v) { + imdsocket * s = (imdsocket *) v; + if (s == NULL) + return; + +#if defined(_MSC_VER) || defined(__MINGW32_VERSION) + closesocket(s->sd); +#else + close(s->sd); +#endif + free(s); +} + +int imdsock_selread(void *v, int sec) { + imdsocket *s = (imdsocket *)v; + fd_set rfd; + struct timeval tv; + int rc; + + if (v == NULL) return 0; + + FD_ZERO(&rfd); + FD_SET(s->sd, &rfd); + memset((void *)&tv, 0, sizeof(struct timeval)); + tv.tv_sec = sec; + do { + rc = select(s->sd+1, &rfd, NULL, NULL, &tv); + } while (rc < 0 && errno == EINTR); + return rc; + +} + +int imdsock_selwrite(void *v, int sec) { + imdsocket *s = (imdsocket *)v; + fd_set wfd; + struct timeval tv; + int rc; + + if (v == NULL) return 0; + + FD_ZERO(&wfd); + FD_SET(s->sd, &wfd); + memset((void *)&tv, 0, sizeof(struct timeval)); + tv.tv_sec = sec; + do { + rc = select(s->sd + 1, NULL, &wfd, NULL, &tv); + } while (rc < 0 && errno == EINTR); + return rc; +} + +/* end of socket code. */ +/*************************************************************************/ + +/*************************************************************************/ +/* start of imd API code. */ +/* Only works with aligned 4-byte quantities, will cause a bus error */ +/* on some platforms if used on unaligned data. */ +void swap4_aligned(void *v, long ndata) { + int *data = (int *) v; + long i; + int *N; + for (i=0; i>24)&0xff) | ((*N&0xff)<<24) | + ((*N>>8)&0xff00) | ((*N&0xff00)<<8)); + } +} + + +/** structure used to perform byte swapping operations */ +typedef union { + int32 i; + struct { + unsigned int highest : 8; + unsigned int high : 8; + unsigned int low : 8; + unsigned int lowest : 8; + } b; +} netint; + +static int32 imd_htonl(int32 h) { + netint n; + n.b.highest = h >> 24; + n.b.high = h >> 16; + n.b.low = h >> 8; + n.b.lowest = 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; + + nleft = n; + while (nleft > 0) { + if ((nread = imdsock_read(s, ptr, nleft)) < 0) { + if (errno == EINTR) + nread = 0; /* and call read() again */ + else + return -1; + } else if (nread == 0) + break; /* EOF */ + nleft -= nread; + ptr += nread; + } + return n-nleft; +} + +static int32 imd_writen(void *s, const char *ptr, int32 n) { + int32 nleft; + int32 nwritten; + + nleft = n; + while (nleft > 0) { + if ((nwritten = imdsock_write(s, ptr, nleft)) <= 0) { + if (errno == EINTR) + nwritten = 0; + else + return -1; + } + nleft -= nwritten; + ptr += nwritten; + } + return n; +} + +int imd_disconnect(void *s) { + IMDheader header; + imd_fill_header(&header, IMD_DISCONNECT, 0); + return (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE); +} + +int imd_handshake(void *s) { + IMDheader header; + imd_fill_header(&header, IMD_HANDSHAKE, 1); + header.length = IMDVERSION; /* Not byteswapped! */ + return (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE); +} + +/* The IMD receive functions */ + +IMDType imd_recv_header(void *s, int32 *length) { + IMDheader header; + if (imd_readn(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE) + return IMD_IOERROR; + swap_header(&header); + *length = header.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); +} + +// Local Variables: +// mode: c++ +// compile-command: "make -j4 openmpi" +// c-basic-offset: 2 +// fill-column: 76 +// indent-tabs-mode: nil +// End: + diff --git a/src/USER-MISC/fix_imd.h b/src/USER-MISC/fix_imd.h new file mode 100644 index 0000000000..a7c3cb36eb --- /dev/null +++ b/src/USER-MISC/fix_imd.h @@ -0,0 +1,131 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + The FixIMD class contains code from VMD and NAMD which is copyrighted + by the Board of Trustees of the University of Illinois and is free to + use with LAMMPS according to point 2 of the UIUC license (10% clause): + +" Licensee may, at its own expense, create and freely distribute +complimentary works that interoperate with the Software, directing others to +the TCBG server to license and obtain the Software itself. Licensee may, at +its own expense, modify the Software to make derivative works. Except as +explicitly provided below, this License shall apply to any derivative work +as it does to the original Software distributed by Illinois. Any derivative +work should be clearly marked and renamed to notify users that it is a +modified version and not the original Software distributed by Illinois. +Licensee agrees to reproduce the copyright notice and other proprietary +markings on any derivative work and to include in the documentation of such +work the acknowledgement: + + "This software includes code developed by the Theoretical and Computational + Biophysics Group in the Beckman Institute for Advanced Science and + Technology at the University of Illinois at Urbana-Champaign." + +Licensee may redistribute without restriction works with up to 1/2 of their +non-comment source code derived from at most 1/10 of the non-comment source +code developed by Illinois and contained in the Software, provided that the +above directions for notice and acknowledgement are observed. Any other +distribution of the Software or any derivative work requires a separate +license with Illinois. Licensee may contact Illinois (vmd@ks.uiuc.edu) to +negotiate an appropriate license for such distribution." +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(imd,FixIMD) + +#else + +#ifndef LMP_FIX_IMD_H +#define LMP_FIX_IMD_H + +#include "fix.h" + +#if defined(LAMMPS_ASYNC_IMD) +#include +#endif + +/* prototype for c wrapper that calls the real worker */ +extern "C" void *fix_imd_ioworker(void *); + +namespace LAMMPS_NS { + +class FixIMD : public Fix { + public: + FixIMD(class LAMMPS *, int, char **); + ~FixIMD(); + int setmask(); + void init(); + void setup(int); + void post_force(int); + void post_force_respa(int, int, int); + double memory_usage(); + + protected: + int imd_port; + void *localsock; + void *clientsock; + + 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 *idmap; // hash for mapping atom indices to consistent order. + int *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_inactive; // true if IMD connection stopped. + int imd_terminate; // true if IMD requests termination of run. + int imd_trate; // IMD transmission rate. + + int unwrap_flag; // true if coordinates need to be unwrapped before sending + 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. + + int me; // my MPI rank in this "world". + int nlevels_respa; // flag to determine respa levels. + + int msglen; + char *msgdata; + +#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 + pthread_cond_t write_cond; // conditional variable for coordinate i/o + pthread_mutex_t read_mutex; // mutex for accessing data receieved by i/o thread + pthread_t iothread; // thread id for i/o thread. + pthread_attr_t iot_attr; // i/o thread attributes. +public: + void ioworker(void); +#endif + +protected: + int reconnect(); +}; + +} + +#endif +#endif + +// Local Variables: +// mode: c++ +// compile-command: "make -j4 openmpi" +// c-basic-offset: 2 +// fill-column: 76 +// indent-tabs-mode: nil +// End: diff --git a/src/USER-MISC/fix_smd.cpp b/src/USER-MISC/fix_smd.cpp new file mode 100644 index 0000000000..4d251d3fe3 --- /dev/null +++ b/src/USER-MISC/fix_smd.cpp @@ -0,0 +1,403 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (UPenn) + based on fix spring by: Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "fix_smd.h" +#include "atom.h" +#include "comm.h" +#include "update.h" +#include "respa.h" +#include "domain.h" +#include "error.h" +#include "group.h" + +using namespace LAMMPS_NS; + +enum { SMD_NONE=0, + SMD_TETHER=1<<0, SMD_COUPLE=1<<1, + SMD_CVEL=1<<2, SMD_CFOR=1<<3, + SMD_AUTOX=1<<4, SMD_AUTOY=1<<5, SMD_AUTOZ=1<<6}; + +#define SMALL 0.001 + +/* ---------------------------------------------------------------------- */ + +FixSMD::FixSMD(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + styleflag = SMD_NONE; + k_smd = f_smd = v_smd = -1.0; + xflag = yflag = zflag = 1; + xc = yc = zc = 0.0; + xn = yn = zn = 1.0; + pmf = r_old = r_now = r0 = 0.0; + + restart_global = 1; + vector_flag = 1; + size_vector = 7; + global_freq = 1; + extvector = 1; + + int argoffs=3; + if (strcmp(arg[argoffs],"cvel") == 0) { + if (narg < argoffs+3) error->all("Illegal fix smd command"); + styleflag |= SMD_CVEL; + k_smd = atof(arg[argoffs+1]); + v_smd = atof(arg[argoffs+2]); // to be multiplied by update->dt when used. + argoffs += 3; + } else if (strcmp(arg[argoffs],"cfor") == 0) { + if (narg < argoffs+2) error->all("Illegal fix smd command"); + styleflag |= SMD_CFOR; + f_smd = atof(arg[argoffs+1]); + argoffs += 2; + } else error->all("Illegal fix smd command"); + + if (strcmp(arg[argoffs],"tether") == 0) { + if (narg < argoffs+5) error->all("Illegal fix smd command"); + styleflag |= SMD_TETHER; + if (strcmp(arg[argoffs+1],"NULL") == 0) xflag = 0; + else xc = atof(arg[argoffs+1]); + if (strcmp(arg[argoffs+2],"NULL") == 0) yflag = 0; + else yc = atof(arg[argoffs+2]); + if (strcmp(arg[argoffs+3],"NULL") == 0) zflag = 0; + else zc = atof(arg[argoffs+3]); + r0 = atof(arg[argoffs+4]); + if (r0 < 0) error->all("R0 < 0 for fix smd command"); + argoffs += 5; + } else if (strcmp(arg[argoffs],"couple") == 0) { + if (narg < argoffs+6) error->all("Illegal fix smd command"); + styleflag |= SMD_COUPLE; + igroup2 = group->find(arg[argoffs+1]); + if (igroup2 == -1) + error->all("Could not find fix smd couple group ID"); + if (igroup2 == igroup) + error->all("Two groups cannot be the same in fix smd couple"); + group2bit = group->bitmask[igroup2]; + + if (strcmp(arg[argoffs+2],"NULL") == 0) xflag = 0; + else if (strcmp(arg[argoffs+2],"auto") == 0) styleflag |= SMD_AUTOX; + else xc = atof(arg[argoffs+2]); + if (strcmp(arg[argoffs+3],"NULL") == 0) yflag = 0; + else if (strcmp(arg[argoffs+3],"auto") == 0) styleflag |= SMD_AUTOY; + else yc = atof(arg[argoffs+3]); + if (strcmp(arg[argoffs+4],"NULL") == 0) zflag = 0; + else if (strcmp(arg[argoffs+4],"auto") == 0) styleflag |= SMD_AUTOZ; + else zc = atof(arg[argoffs+4]); + + r0 = atof(arg[argoffs+5]); + if (r0 < 0) error->all("R0 < 0 for fix smd command"); + argoffs +=6; + } else error->all("Illegal fix smd command"); + + force_flag = 0; + ftotal[0] = ftotal[1] = ftotal[2] = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +int FixSMD::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= POST_FORCE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::init() +{ + double xcm[3], xcm2[3]; + masstotal = group->mass(igroup); + group->xcm(igroup,masstotal,xcm); + + double dx,dy,dz; + if (styleflag & SMD_TETHER) { + dx = xc - xcm[0]; + dy = yc - xcm[1]; + dz = zc - xcm[2]; + } else { /* SMD_COUPLE */ + masstotal2 = group->mass(igroup2); + group->xcm(igroup2,masstotal2,xcm2); + if (styleflag & SMD_AUTOX) dx = xcm2[0] - xcm[0]; + else dx = xc; + if (styleflag & SMD_AUTOY) dy = xcm2[1] - xcm[1]; + else dy = yc; + if (styleflag & SMD_AUTOZ) dz = xcm2[2] - xcm[2]; + else dz = zc; + } + + if (!xflag) dx = 0.0; + if (!yflag) dy = 0.0; + if (!zflag) dz = 0.0; + r_old = sqrt(dx*dx + dy*dy + dz*dz); + if (r_old > SMALL) { + xn = dx/r_old; + yn = dy/r_old; + zn = dz/r_old; + } + + if (strstr(update->integrate_style,"respa")) + nlevels_respa = ((Respa *) update->integrate)->nlevels; +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::setup(int vflag) +{ + if (strstr(update->integrate_style,"verlet")) + post_force(vflag); + else { + ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1); + post_force_respa(vflag,nlevels_respa-1,0); + ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::post_force(int vflag) +{ + if (styleflag & SMD_TETHER) smd_tether(); + else smd_couple(); + + if (styleflag & SMD_CVEL) r_old += v_smd * update->dt; +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::smd_tether() +{ + double xcm[3]; + group->xcm(igroup,masstotal,xcm); + + // fx,fy,fz = components of k * (r-r0) + + double dx,dy,dz,fx,fy,fz,r,dr; + + dx = xcm[0] - xc; + dy = xcm[1] - yc; + dz = xcm[2] - zc; + r_now = sqrt(dx*dx + dy*dy + dz*dz); + + if (!xflag) dx = 0.0; + if (!yflag) dy = 0.0; + if (!zflag) dz = 0.0; + r = sqrt(dx*dx + dy*dy + dz*dz); + if (styleflag & SMD_CVEL) { + if(r > SMALL) { + double fsign; + fsign = (v_smd<0.0) ? -1.0 : 1.0; + dr = r - r0 - r_old; + fx = k_smd*dx*dr/r; + fy = k_smd*dy*dr/r; + fz = k_smd*dz*dr/r; + pmf += (fx*xn + fy*yn + fz*zn) * v_smd * update->dt; + } + } else { + r_old = r; + fx = f_smd*dx/r; + fy = f_smd*dy/r; + fz = f_smd*dz/r; + } + + // apply restoring force to atoms in group + // f = -k*(r-r0)*mass/masstotal + + double **f = atom->f; + int *mask = atom->mask; + int *type = atom->type; + double *mass = atom->mass; + int nlocal = atom->nlocal; + + ftotal[0] = ftotal[1] = ftotal[2] = 0.0; + force_flag = 0; + + double massfrac; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + massfrac = mass[type[i]]/masstotal; + f[i][0] -= fx*massfrac; + f[i][1] -= fy*massfrac; + f[i][2] -= fz*massfrac; + ftotal[0] -= fx*massfrac; + ftotal[1] -= fy*massfrac; + ftotal[2] -= fz*massfrac; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::smd_couple() +{ + double xcm[3],xcm2[3]; + group->xcm(igroup,masstotal,xcm); + group->xcm(igroup2,masstotal2,xcm2); + + // renormalize direction of spring + double dx,dy,dz,r,dr; + if (styleflag & SMD_AUTOX) dx = xcm2[0] - xcm[0]; + else dx = xn*r_old; + if (styleflag & SMD_AUTOY) dy = xcm2[1] - xcm[1]; + else dy = yn*r_old; + if (styleflag & SMD_AUTOZ) dz = xcm2[2] - xcm[2]; + else dz = zn*r_old; + if (!xflag) dx = 0.0; + if (!yflag) dy = 0.0; + if (!zflag) dz = 0.0; + r = sqrt(dx*dx + dy*dy + dz*dz); + if (r > SMALL) { + xn = dx/r; yn = dy/r; zn = dz/r; + } + + double fx,fy,fz; + if (styleflag & SMD_CVEL) { + dx = xcm2[0] - xcm[0]; + dy = xcm2[1] - xcm[1]; + dz = xcm2[2] - xcm[2]; + r_now = sqrt(dx*dx + dy*dy + dz*dz); + + dx -= xn*r_old; + dy -= yn*r_old; + dz -= zn*r_old; + + if (!xflag) dx = 0.0; + if (!yflag) dy = 0.0; + if (!zflag) dz = 0.0; + r = sqrt(dx*dx + dy*dy + dz*dz); + dr = r - r0; + + if (r > SMALL) { + double fsign; + fsign = (v_smd<0.0) ? -1.0 : 1.0; + + fx = k_smd*dx*dr/r; + fy = k_smd*dy*dr/r; + fz = k_smd*dz*dr/r; + pmf += (fx*xn + fy*yn + fz*zn) * fsign * v_smd * update->dt; + } + + } else { + dx = xcm2[0] - xcm[0]; + dy = xcm2[1] - xcm[1]; + dz = xcm2[2] - xcm[2]; + r_now = sqrt(dx*dx + dy*dy + dz*dz); + r_old = r; + + fx = f_smd*xn; + fy = f_smd*yn; + fz = f_smd*zn; + } + + // apply restoring force to atoms in group + // f = -k*(r-r0)*mass/masstotal + + double **f = atom->f; + int *mask = atom->mask; + int *type = atom->type; + double *mass = atom->mass; + int nlocal = atom->nlocal; + + ftotal[0] = ftotal[1] = ftotal[2] = 0.0; + force_flag = 0; + + double massfrac; + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + massfrac = mass[type[i]]/masstotal; + f[i][0] += fx*massfrac; + f[i][1] += fy*massfrac; + f[i][2] += fz*massfrac; + ftotal[0] += fx*massfrac; + ftotal[1] += fy*massfrac; + ftotal[2] += fz*massfrac; + } + if (mask[i] & group2bit) { + massfrac = mass[type[i]]/masstotal2; + f[i][0] -= fx*massfrac; + f[i][1] -= fy*massfrac; + f[i][2] -= fz*massfrac; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::write_restart(FILE *fp) +{ +#define RESTART_ITEMS 5 + double buf[RESTART_ITEMS], fsign; + + if (comm->me == 0) { + // make sure we project the force into the direction of the pulling. + fsign = (v_smd<0.0) ? -1.0 : 1.0; + buf[0] = r_old; + buf[1] = xn*fsign; + buf[2] = yn*fsign; + buf[3] = zn*fsign; + buf[4] = pmf; + int size = RESTART_ITEMS*sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(&buf[0],sizeof(double),RESTART_ITEMS,fp); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::restart(char *buf) +{ + double *list = (double *)buf; + r_old = list[0]; + xn=list[1]; + yn=list[2]; + zn=list[3]; + pmf=list[4]; +} + +/* ---------------------------------------------------------------------- */ + +void FixSMD::post_force_respa(int vflag, int ilevel, int iloop) +{ + if (ilevel == nlevels_respa-1) post_force(vflag); +} + +/* ---------------------------------------------------------------------- + return components of total smd force on fix group +------------------------------------------------------------------------- */ + +double FixSMD::compute_vector(int n) +{ + // only sum across procs one time + + if (force_flag == 0) { + MPI_Allreduce(ftotal,ftotal_all,3,MPI_DOUBLE,MPI_SUM,world); + force_flag = 1; + if (styleflag & SMD_CVEL) { + ftotal_all[3]=ftotal_all[0]*xn+ftotal_all[1]*yn+ftotal_all[2]*zn; + ftotal_all[4]=r_old; + } else { + ftotal_all[3]=f_smd; + ftotal_all[4]=r_old; + } + ftotal_all[5]=r_now; + ftotal_all[6]=pmf; + } + return ftotal_all[n]; +} diff --git a/src/USER-MISC/fix_smd.h b/src/USER-MISC/fix_smd.h new file mode 100644 index 0000000000..5319ca9558 --- /dev/null +++ b/src/USER-MISC/fix_smd.h @@ -0,0 +1,60 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(smd,FixSMD) + +#else + +#ifndef LMP_FIX_SMD_H +#define LMP_FIX_SMD_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixSMD : public Fix { + public: + FixSMD(class LAMMPS *, int, char **); + int setmask(); + void init(); + void setup(int); + void post_force(int); + void post_force_respa(int, int, int); + double compute_vector(int); + + void write_restart(FILE *); + void restart(char *); + + private: + double xc,yc,zc,xn,yn,zn,r0; + double k_smd,f_smd,v_smd; + int xflag,yflag,zflag; + int styleflag; + double r_old,r_now,pmf; + + int igroup2,group2bit; + double masstotal,masstotal2; + int nlevels_respa; + double ftotal[3],ftotal_all[7]; + int force_flag; + + void smd_tether(); + void smd_couple(); +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/pair_cdeam.cpp b/src/USER-MISC/pair_cdeam.cpp new file mode 100644 index 0000000000..2d0e224f1c --- /dev/null +++ b/src/USER-MISC/pair_cdeam.cpp @@ -0,0 +1,642 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Alexander Stukowski + Technical University of Darmstadt, + Germany Department of Materials Science +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_cdeam.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +// This is for debugging purposes. The ASSERT() macro is used in the code to check +// if everything runs as expected. Change this to #if 0 if you don't need the checking. +#if 0 + #define ASSERT(cond) ((!(cond)) ? my_failure(error,__FILE__,__LINE__) : my_noop()) + + inline void my_noop() {} + inline void my_failure(Error* error, const char* file, int line) { + char str[1024]; + sprintf(str,"Assertion failure: File %s, line %i", file, line); + error->one(str); + } +#else + #define ASSERT(cond) +#endif + +#define MAXLINE 1024 // This sets the maximum line length in EAM input files. + +PairCDEAM::PairCDEAM(LAMMPS *lmp, int _cdeamVersion) : PairEAM(lmp), PairEAMAlloy(lmp), cdeamVersion(_cdeamVersion) +{ + single_enable = 0; + rhoB = NULL; + D_values = NULL; + hcoeff = NULL; + + // Set communication buffer sizes needed by this pair style. + if(cdeamVersion == 1) { + comm_forward = 4; + comm_reverse = 3; + } + else if(cdeamVersion == 2) { + comm_forward = 3; + comm_reverse = 2; + } + else { + error->all("Invalid CD-EAM potential version."); + } +} + +PairCDEAM::~PairCDEAM() +{ + memory->destroy(rhoB); + memory->destroy(D_values); + if(hcoeff) delete[] hcoeff; +} + +void PairCDEAM::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,rhoip,rhojp,recip,phi; + int *ilist,*jlist,*numneigh,**firstneigh; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = eflag_global = eflag_atom = 0; + + // Grow per-atom arrays if necessary + if(atom->nmax > nmax) { + memory->destroy(rho); + memory->destroy(fp); + memory->destroy(rhoB); + memory->destroy(D_values); + nmax = atom->nmax; + memory->create(rho,nmax,"pair:rho"); + memory->create(rhoB,nmax,"pair:rhoB"); + memory->create(fp,nmax,"pair:fp"); + memory->create(D_values,nmax,"pair:D_values"); + } + + double **x = atom->x; + double **f = atom->f; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // Zero out per-atom arrays. + int m = nlocal + atom->nghost; + for(i = 0; i < m; i++) { + rho[i] = 0.0; + rhoB[i] = 0.0; + D_values[i] = 0.0; + } + + // Stage I + + // Compute rho and rhoB at each local atom site. + // Additionally calculate the D_i values here if we are using the one-site formulation. + // For the two-site formulation we have to calculate the D values in an extra loop (Stage II). + for(ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for(jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if(rsq < cutforcesq) { + jtype = type[j]; + double r = sqrt(rsq); + const EAMTableIndex index = radiusToTableIndex(r); + double localrho = RhoOfR(index, jtype, itype); + rho[i] += localrho; + if(jtype == speciesB) rhoB[i] += localrho; + if(newton_pair || j < nlocal) { + localrho = RhoOfR(index, itype, jtype); + rho[j] += localrho; + if(itype == speciesB) rhoB[j] += localrho; + } + + if(cdeamVersion == 1 && itype != jtype) { + // Note: if the i-j interaction is not concentration dependent (because either + // i or j are not species A or B) then its contribution to D_i and D_j should + // be ignored. + // This if-clause is only required for a ternary. + if((itype == speciesA && jtype == speciesB) || (jtype == speciesA && itype == speciesB)) { + double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r); + D_values[i] += Phi_AB; + if(newton_pair || j < nlocal) + D_values[j] += Phi_AB; + } + } + } + } + } + + // Communicate and sum densities. + if(newton_pair) { + communicationStage = 1; + comm->reverse_comm_pair(this); + } + + // fp = derivative of embedding energy at each atom + // phi = embedding energy at each atom + for(ii = 0; ii < inum; ii++) { + i = ilist[ii]; + EAMTableIndex index = rhoToTableIndex(rho[i]); + fp[i] = FPrimeOfRho(index, type[i]); + if(eflag) { + phi = FofRho(index, type[i]); + if (eflag_global) eng_vdwl += phi; + if (eflag_atom) eatom[i] += phi; + } + } + + // Communicate derivative of embedding function and densities + // and D_values (this for one-site formulation only). + communicationStage = 2; + comm->forward_comm_pair(this); + + // The electron densities may not drop to zero because then the concentration would no longer be defined. + // But the concentration is not needed anyway if there is no interaction with another atom, which is the case + // if the electron density is exactly zero. That's why the following lines have been commented out. + // + //for(i = 0; i < nlocal + atom->nghost; i++) { + // if(rho[i] == 0 && (type[i] == speciesA || type[i] == speciesB)) + // error->one("CD-EAM potential routine: Detected atom with zero electron density."); + //} + + // Stage II + // This is only required for the original two-site formulation of the CD-EAM potential. + + if(cdeamVersion == 2) { + // Compute intermediate value D_i for each atom. + for(ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // This code line is required for ternary alloys. + if(itype != speciesA && itype != speciesB) continue; + + double x_i = rhoB[i] / rho[i]; // Concentration at atom i. + + for(jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtype = type[j]; + if(itype == jtype) continue; + + // This code line is required for ternary alloys. + if(jtype != speciesA && jtype != speciesB) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if(rsq < cutforcesq) { + double r = sqrt(rsq); + const EAMTableIndex index = radiusToTableIndex(r); + + // The concentration independent part of the cross pair potential. + double Phi_AB = PhiOfR(index, itype, jtype, 1.0 / r); + + // Average concentration of two sites + double x_ij = 0.5 * (x_i + rhoB[j]/rho[j]); + + // Calculate derivative of h(x_ij) polynomial function. + double h_prime = evalHprime(x_ij); + + D_values[i] += h_prime * Phi_AB / (2.0 * rho[i] * rho[i]); + if(newton_pair || j < nlocal) + D_values[j] += h_prime * Phi_AB / (2.0 * rho[j] * rho[j]); + } + } + } + + // Communicate and sum D values. + if(newton_pair) { + communicationStage = 3; + comm->reverse_comm_pair(this); + } + communicationStage = 4; + comm->forward_comm_pair(this); + } + + // Stage III + + // Compute force acting on each atom. + for(ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // Concentration at site i + double x_i = -1.0; // The value -1 indicates: no concentration dependence for all interactions of atom i. + // It will be replaced by the concentration at site i if atom i is either A or B. + + double D_i, h_prime_i; + + // This if-clause is only required for ternary alloys. + if((itype == speciesA || itype == speciesB) && rho[i] != 0.0) { + + // Compute local concentration at site i. + x_i = rhoB[i]/rho[i]; + ASSERT(x_i >= 0 && x_i<=1.0); + + if(cdeamVersion == 1) { + // Calculate derivative of h(x_i) polynomial function. + h_prime_i = evalHprime(x_i); + D_i = D_values[i] * h_prime_i / (2.0 * rho[i] * rho[i]); + } + else if(cdeamVersion == 2) { + D_i = D_values[i]; + } + else ASSERT(false); + } + + for(jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + + if(rsq < cutforcesq) { + jtype = type[j]; + double r = sqrt(rsq); + const EAMTableIndex index = radiusToTableIndex(r); + + // rhoip = derivative of (density at atom j due to atom i) + // rhojp = derivative of (density at atom i due to atom j) + // psip needs both fp[i] and fp[j] terms since r_ij appears in two + // terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji) + // hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip + rhoip = RhoPrimeOfR(index, itype, jtype); + rhojp = RhoPrimeOfR(index, jtype, itype); + fpair = fp[i]*rhojp + fp[j]*rhoip; + recip = 1.0/r; + + double x_j = -1; // The value -1 indicates: no concentration dependence for this i-j pair + // because atom j is not of species A nor B. + + // This code line is required for ternary alloy. + if(jtype == speciesA || jtype == speciesB) { + ASSERT(rho[i] != 0.0); + ASSERT(rho[j] != 0.0); + + // Compute local concentration at site j. + x_j = rhoB[j]/rho[j]; + ASSERT(x_j >= 0 && x_j<=1.0); + + double D_j; + if(cdeamVersion == 1) { + // Calculate derivative of h(x_j) polynomial function. + double h_prime_j = evalHprime(x_j); + D_j = D_values[j] * h_prime_j / (2.0 * rho[j] * rho[j]); + } + else if(cdeamVersion == 2) { + D_j = D_values[j]; + } + else ASSERT(false); + + double t2 = -rhoB[j]; + if(itype == speciesB) t2 += rho[j]; + fpair += D_j * rhoip * t2; + } + + // This if-clause is only required for a ternary alloy. + // Actually we don't need it at all because D_i should be zero anyway if + // atom i has no concentration dependent interactions (because it is not species A or B). + if(x_i != -1.0) { + double t1 = -rhoB[i]; + if(jtype == speciesB) t1 += rho[i]; + fpair += D_i * rhojp * t1; + } + + double phip; + double phi = PhiOfR(index, itype, jtype, recip, phip); + if(itype == jtype || x_i == -1.0 || x_j == -1.0) { + // Case of no concentration dependence. + fpair += phip; + } + else { + // We have a concentration dependence for the i-j interaction. + double h; + if(cdeamVersion == 1) { + // Calculate h(x_i) polynomial function. + double h_i = evalH(x_i); + // Calculate h(x_j) polynomial function. + double h_j = evalH(x_j); + h = 0.5 * (h_i + h_j); + } + else if(cdeamVersion == 2) { + // Average concentration. + double x_ij = 0.5 * (x_i + x_j); + // Calculate h(x_ij) polynomial function. + h = evalH(x_ij); + } + else ASSERT(false); + fpair += h * phip; + phi *= h; + } + + // Divide by r_ij and negate to get forces from gradient. + fpair /= -r; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if(newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if(eflag) evdwl = phi; + if(evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz); + } + } + } + + if(vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairCDEAM::coeff(int narg, char **arg) +{ + PairEAMAlloy::coeff(narg, arg); + + // Make sure the EAM file is a CD-EAM binary alloy. + if(setfl->nelements < 2) + error->all("The EAM file must contain at least 2 elements to be used with the eam/cd pair style."); + + // Read in the coefficients of the h polynomial from the end of the EAM file. + read_h_coeff(arg[2]); + + // Determine which atom type is the A species and which is the B species in the alloy. + // By default take the first element (index 0) in the EAM file as the A species + // and the second element (index 1) in the EAM file as the B species. + speciesA = -1; + speciesB = -1; + for(int i = 1; i <= atom->ntypes; i++) { + if(map[i] == 0) { + if(speciesA >= 0) + error->all("The first element from the EAM file may only be mapped to a single atom type."); + speciesA = i; + } + if(map[i] == 1) { + if(speciesB >= 0) + error->all("The second element from the EAM file may only be mapped to a single atom type."); + speciesB = i; + } + } + if(speciesA < 0) + error->all("The first element from the EAM file must be mapped to exactly one atom type."); + if(speciesB < 0) + error->all("The second element from the EAM file must be mapped to exactly one atom type."); +} + +/* ---------------------------------------------------------------------- + Reads in the h(x) polynomial coefficients +------------------------------------------------------------------------- */ +void PairCDEAM::read_h_coeff(char *filename) +{ + if(comm->me == 0) { + // Open potential file + FILE *fp; + char line[MAXLINE]; + char nextline[MAXLINE]; + fp = fopen(filename,"r"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open EAM potential file %s", filename); + error->one(str); + } + + // h coefficients are stored at the end of the file. + // Skip to last line of file. + while(fgets(nextline, MAXLINE, fp) != NULL) { + strcpy(line, nextline); + } + char* ptr = strtok(line, " \t\n\r\f"); + int degree = atoi(ptr); + nhcoeff = degree+1; + hcoeff = new double[nhcoeff]; + int i = 0; + while((ptr = strtok(NULL," \t\n\r\f")) != NULL && i < nhcoeff) { + hcoeff[i++] = atof(ptr); + } + if(i != nhcoeff || nhcoeff < 1) + error->one("Failed to read h(x) function coefficients from EAM file."); + + // Close the potential file. + fclose(fp); + } + + MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world); + if(comm->me != 0) hcoeff = new double[nhcoeff]; + MPI_Bcast(hcoeff, nhcoeff, MPI_DOUBLE, 0, world); +} + + +/* ---------------------------------------------------------------------- */ + +int PairCDEAM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) +{ + int i,j,m; + + m = 0; + if(communicationStage == 2) { + if(cdeamVersion == 1) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = fp[j]; + buf[m++] = rho[j]; + buf[m++] = rhoB[j]; + buf[m++] = D_values[j]; + } + return 4; + } + else if(cdeamVersion == 2) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = fp[j]; + buf[m++] = rho[j]; + buf[m++] = rhoB[j]; + } + return 3; + } + else { ASSERT(false); return 0; } + } + else if(communicationStage == 4) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = D_values[j]; + } + return 1; + } + else return 0; +} + +/* ---------------------------------------------------------------------- */ + +void PairCDEAM::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + if(communicationStage == 2) { + if(cdeamVersion == 1) { + for(i = first; i < last; i++) { + fp[i] = buf[m++]; + rho[i] = buf[m++]; + rhoB[i] = buf[m++]; + D_values[i] = buf[m++]; + } + } + else if(cdeamVersion == 2) { + for(i = first; i < last; i++) { + fp[i] = buf[m++]; + rho[i] = buf[m++]; + rhoB[i] = buf[m++]; + } + } + else ASSERT(false); + } + else if(communicationStage == 4) { + for(i = first; i < last; i++) { + D_values[i] = buf[m++]; + } + } +} + +/* ---------------------------------------------------------------------- */ +int PairCDEAM::pack_reverse_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + + if(communicationStage == 1) { + if(cdeamVersion == 1) { + for(i = first; i < last; i++) { + buf[m++] = rho[i]; + buf[m++] = rhoB[i]; + buf[m++] = D_values[i]; + } + return 3; + } + else if(cdeamVersion == 2) { + for(i = first; i < last; i++) { + buf[m++] = rho[i]; + buf[m++] = rhoB[i]; + } + return 2; + } + else { ASSERT(false); return 0; } + } + else if(communicationStage == 3) { + for(i = first; i < last; i++) { + buf[m++] = D_values[i]; + } + return 1; + } + else return 0; +} + +/* ---------------------------------------------------------------------- */ + +void PairCDEAM::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,j,m; + + m = 0; + if(communicationStage == 1) { + if(cdeamVersion == 1) { + for(i = 0; i < n; i++) { + j = list[i]; + rho[j] += buf[m++]; + rhoB[j] += buf[m++]; + D_values[j] += buf[m++]; + } + } + else if(cdeamVersion == 2) { + for(i = 0; i < n; i++) { + j = list[i]; + rho[j] += buf[m++]; + rhoB[j] += buf[m++]; + } + } + else ASSERT(false); + } + else if(communicationStage == 3) { + for(i = 0; i < n; i++) { + j = list[i]; + D_values[j] += buf[m++]; + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ +double PairCDEAM::memory_usage() +{ + double bytes = 2 * nmax * sizeof(double); + return PairEAMAlloy::memory_usage() + bytes; +} diff --git a/src/USER-MISC/pair_cdeam.h b/src/USER-MISC/pair_cdeam.h new file mode 100644 index 0000000000..e39923d568 --- /dev/null +++ b/src/USER-MISC/pair_cdeam.h @@ -0,0 +1,230 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(eam/cd,PairCDEAM_OneSite) +PairStyle(eam/cd/old,PairCDEAM_TwoSite) + +#else + +#ifndef LMP_PAIR_CDEAM_H +#define LMP_PAIR_CDEAM_H + +#include "pair_eam_alloy.h" + +namespace LAMMPS_NS { + +class PairCDEAM : public PairEAMAlloy +{ +public: + /// Constructor. + PairCDEAM(class LAMMPS*, int cdeamVersion); + + /// Destructor. + virtual ~PairCDEAM(); + + /// Calculates the energies and forces for all atoms in the system. + void compute(int, int); + + /// Parses the pair_coeff command parameters for this pair style. + void coeff(int, char **); + + /// This is for MPI communication with neighbor nodes. + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + int pack_reverse_comm(int, int, double *); + void unpack_reverse_comm(int, int *, double *); + + /// Reports the memory usage of this pair style to LAMMPS. + double memory_usage(); + + /// Parses the coefficients of the h polynomial from the end of the EAM file. + void read_h_coeff(char* filename); + + public: + // The public interface exposed by this potential class. + + // Evaluates the h(x) polynomial for a given local concentration x. + inline double evalH(double x) const { + double v = 0.0; + for(int i = nhcoeff-1; i >= 1; i--) { + v = (v + hcoeff[i]) * x; + } + return v + hcoeff[0]; + } + + // Calculates the derivative of the h(x) polynomial. + inline double evalHprime(double x) const { + double v = 0.0; + for(int i = nhcoeff-1; i >= 2; i--) { + v = (v + (double)i * hcoeff[i]) * x; + } + return v + hcoeff[1]; + } + + // We have two versions of the CD-EAM potential. The user can choose which one he wants to use: + // + // Version 1 - One-site concentration: The h(x_i) function depends only on the concentration at the atomic site i. + // This is a new version with a slight modification to the formula. It happens to be computationally more efficient. + // It has been published in the MSMSE 2009 paper. + // + // Version 2 - Two-site concentration: The h(x_ij) function depends on the concentrations at two atomic sites i and j. + // This is the original version from the 2005 PRL paper. + int cdeamVersion; + + // Per-atom arrays + + // The partial density of B atoms at each atom site. + double *rhoB; + + // The intermediate value D_i for each atom. + // The meaning of these values depend on the version of the CD-EAM potential used: + // + // For the one-site concentration version these are the v_i values defined in equation (21) + // of the MSMSE paper. + // + // For the old two-site concentration version these are the M_i values defined in equation (12) + // of the MSMSE paper. + double *D_values; + + // The atom type index that is considered to be the A species in the alloy. + int speciesA; + + // The atom type index that is considered to be the B species in the alloy. + int speciesB; + + protected: + + // Evaluation functions: + + // This structure specifies an entry in one of the EAM spline tables + // and the corresponding floating point part. + typedef struct { + int m; + double p; + } EAMTableIndex; + + // Converts a radius value to an index value to be used in a spline table lookup. + inline EAMTableIndex radiusToTableIndex(double r) const { + EAMTableIndex index; + index.p = r*rdr + 1.0; + index.m = static_cast(index.p); + index.m = index.m <= (nr-1) ? index.m : (nr-1); + index.p -= index.m; + index.p = index.p <= 1.0 ? index.p : 1.0; + return index; + } + + // Converts a density value to an index value to be used in a spline table lookup. + inline EAMTableIndex rhoToTableIndex(double rho) const { + EAMTableIndex index; + index.p = rho*rdrho + 1.0; + index.m = static_cast(index.p); + index.m = index.m <= (nrho-1) ? index.m : (nrho-1); + index.p -= index.m; + index.p = index.p <= 1.0 ? index.p : 1.0; + return index; + } + + // Computes the derivative of rho(r) + inline double RhoPrimeOfR(const EAMTableIndex& index, int itype, int jtype) const { + const double* coeff = rhor_spline[type2rhor[itype][jtype]][index.m]; + return (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; + } + + // Computes rho(r) + inline double RhoOfR(const EAMTableIndex& index, int itype, int jtype) const { + const double* coeff = rhor_spline[type2rhor[itype][jtype]][index.m]; + return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; + } + + // Computes the derivative of F(rho) + inline double FPrimeOfRho(const EAMTableIndex& index, int itype) const { + const double* coeff = frho_spline[type2frho[itype]][index.m]; + return (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; + } + + // Computes F(rho) + inline double FofRho(const EAMTableIndex& index, int itype) const { + const double* coeff = frho_spline[type2frho[itype]][index.m]; + return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; + } + + // Computes the derivative of z2(r) + inline double Z2PrimeOfR(const EAMTableIndex& index, int itype, int jtype) const { + const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; + return (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; + } + + // Computes z2(r) + inline double Z2OfR(const EAMTableIndex& index, int itype, int jtype) const { + const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; + return ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; + } + + // Computes pair potential V_ij(r). + inline double PhiOfR(const EAMTableIndex& index, int itype, int jtype, const double oneOverR) const { + // phi = pair potential energy + // z2 = phi * r + const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; + const double z2 = ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; + return z2 * oneOverR; + } + + // Computes pair potential V_ij(r) and its derivative. + inline double PhiOfR(const EAMTableIndex& index, int itype, int jtype, const double oneOverR, double& phid) const { + // phi = pair potential energy + // phip = phi' + // z2 = phi * r + // z2p = (phi * r)' = (phi' r) + phi + const double* coeff = z2r_spline[type2z2r[itype][jtype]][index.m]; + const double z2p = (coeff[0]*index.p + coeff[1])*index.p + coeff[2]; + const double z2 = ((coeff[3]*index.p + coeff[4])*index.p + coeff[5])*index.p + coeff[6]; + const double phi = z2 * oneOverR; + phid = z2p * oneOverR - phi * oneOverR; + return phi; + } + + // Parameters + + // h() polynomial function coefficients + double* hcoeff; + // The number of coefficients in the polynomial. + int nhcoeff; + + // This specifies the calculation stage to let the pack/unpack communication routines know + // which data to exchange. + int communicationStage; +}; + +/// The one-site concentration formulation of CD-EAM. + class PairCDEAM_OneSite : public PairCDEAM + { + public: + /// Constructor. + PairCDEAM_OneSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 1) {} + }; + + /// The two-site concentration formulation of CD-EAM. + class PairCDEAM_TwoSite : public PairCDEAM + { + public: + /// Constructor. + PairCDEAM_TwoSite(class LAMMPS* lmp) : PairEAM(lmp), PairCDEAM(lmp, 2) {} + }; + +}; + +#endif +#endif