From: Dalhaimer, Paul Matthew (pdalhaim_at_utk.edu)
Date: Fri Apr 30 2021 - 06:24:39 CDT
Another option, for those who know C++, is to manually change the source code in GlobalMasterSMD.C:
/**
*** Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by
*** The Board of Trustees of the University of Illinois.
*** All rights reserved.
**/
#include "InfoStream.h"
#include "Node.h"
#include "Molecule.h"
#include "NamdTypes.h"
#include "GlobalMaster.h"
#include "GlobalMasterSMD.h"
#include "PDB.h"
#include "PDBData.h"
/* XXX necessary?
#include <iostream.h>
#if !defined(WIN32) || defined(__CYGWIN__)
#include <strstream.h>
#else
#include <strstrea.h>
#endif
#include <string.h>
*/
//#define DEBUGM
#define MIN_DEBUG_LEVEL 1
#include "Debug.h"
GlobalMasterSMD::GlobalMasterSMD(BigReal spring_constant,
BigReal transverse_spring_constant, BigReal velocity,
const Vector direction, int output_frequency,
int first_timestep, const char *filename,
int numAtoms) {
DebugM(3,"initialize called\n");
moveVel = velocity;
moveDir = direction;
outputFreq = output_frequency;
k = spring_constant;
k2 = transverse_spring_constant;
currentTime = first_timestep;
parseAtoms(filename,numAtoms);
iout << iINFO << requestedGroups()[0].size() << " SMD ATOMS\n" << endi;
DebugM(1,"done with initialize\n");
}
void GlobalMasterSMD::parseAtoms(const char *file, int numTotalAtoms) {
DebugM(3,"parseAtoms called\n");
PDB smdpdb(file);
int numatoms = smdpdb.num_atoms();
if (numatoms < 1)
NAMD_die("No atoms found in SMDFile\n");
if (numatoms != numTotalAtoms)
NAMD_die("The number of atoms in SMDFile must be equal to the total number of atoms in the structure!");
// add a single group
modifyRequestedGroups().resize(1);
modifyGroupForces().resize(1);
// add all the required atoms to the group
// also compute the center of mass of these atoms
cm.x = cm.y = cm.z = 0;
Molecule *mol = Node::Object()->molecule; // to get masses
BigReal imass = 0;
for (int i=0; i<numatoms; i++) {
#ifdef MEM_OPT_VERSION
PDBCoreData *atom = smdpdb.atom(i);
#else
PDBAtom *atom = smdpdb.atom(i); // get an atom from the file
#endif
if (atom->occupancy()) { // if occupancy is not 0, then add it!
// add the atom to the list
modifyRequestedGroups()[0].add(i);
// compute the center of mass
BigReal mass = mol->atommass(i);
cm.x += atom->xcoor()*mass;
cm.y += atom->ycoor()*mass;
cm.z += atom->zcoor()*mass;
imass += mass;
}
}
if (imass == 0) // we didn't find any!
NAMD_die("SMDFile contained no SMD atoms (atoms w/ nonzero occupancy)\n");
// now compute the center of mass
cm /= imass;
DebugM(1,"done with parseAtoms\n");
}
GlobalMasterSMD::~GlobalMasterSMD() { }
void GlobalMasterSMD::calculate() {
DebugM(3,"calculate called\n");
// make extra-sure we have been set up correctly
if(!requestedGroups().size() == 1)
NAMD_die("Internal error: uninitialized!");
// set the correct forces
Position curcm = *getGroupPositionBegin(); // get the center of mass
DebugM(1,"Current CM "<<cm<<"\n");
//BigReal diff = (curcm - cm)*moveDir;
// second term below is along transverse direction: -diff*moveDir + (diff*moveDir - (curcm-cm)) = -(curcm-cm)
// so if k = k2 and moveVel = 0, we see that f = -k * (curcm - cm), the desired result
// Force f = k*(moveVel*currentTime - diff)*moveDir + k2*(diff*moveDir - (curcm - cm));
Force f = k*(moveDir-cumcm);
modifyGroupForces()[0] = f;
// print some output sometimes
if (currentTime % outputFreq == 0)
output(currentTime, curcm, f);
// keep track of the number of timesteps
currentTime++;
DebugM(1,"done with calculate: force: "<<f<<"\n");
}
void GlobalMasterSMD::output(int t, Position p, Force f) {
DebugM(3,"output called\n");
if (t % (100*outputFreq) == 0)
iout << "SMDTITLE: TS CURRENT_POSITION FORCE\n" << endi;
iout << "SMD CBE 579 " << t << ' ' << p << ' ' << f*PNPERKCALMOL << '\n' << endi;
DebugM(1,"done with output\n");
}
From: owner-namd-l_at_ks.uiuc.edu <owner-namd-l_at_ks.uiuc.edu> on behalf of Hemanth Haridas <hemanth.h_at_iitgn.ac.in>
Date: Friday, April 30, 2021 at 7:17 AM
To: namd-l_at_ks.uiuc.edu <namd-l_at_ks.uiuc.edu>, SHIVAM TIWARI <t.shivam_at_iitg.ac.in>
Subject: Re: namd-l: graphene sheet as a rigid piston
Shivam,
I had some luck with using the SMD module to get graphene working as a piston. The trick is to tag one graphene layer as moving (SMD active), and to restrain the layer containing the nanopore with a very high harmonic restraint.
Hope this helps,
Regards
Hemanth
On Wed, Apr 28, 2021 at 11:13 PM SHIVAM TIWARI <t.shivam_at_iitg.ac.in<mailto:t.shivam_at_iitg.ac.in>> wrote:
Dear all,
I want to simulate reverse osmosis process on a box of water and ions across a membrane, for which, I am applying a transmembrane pressure on the water box through a graphene sheet (as a piston). I am using tclforces to achieve this, the script I am using is as follows:
tclForces on
tclForcesScript {
set linaccel1 {0.0 0.0 -0.0006449914635}
set linaccel2 {0.0 0.0 0.0006449914635}
set watIdList1 {}
set watIdList2 {}
for {set i 1} {$i <= 5520} {incr i} {
lappend watIdList1 $i
addatom $i
}
for {set i 229378} {$i <= 234897} {incr i} {
lappend watIdList2 $i
addatom $i
}
But the problem with this is that the graphene sheet behaves as a graphene sheet and not as a rigid piston, that is, the graphene sheet gets deformed (curled up) as soon as the pressure is applied. Although, the water does get pushed, but I suspect that the pressure applied in this case is much lesser than the desired, since most of the force applied on the graphene sheet must be getting wasted on the deformation of the graphene sheet. So how can I achieve a rigid graphene piston?
regards
shivam
-- Hemanth H 18310019 Research Scholar IIT Gandhinagar
This archive was generated by hypermail 2.1.6 : Fri Dec 31 2021 - 23:17:11 CST