Skip to first unread message
Marc Riera
Jun 2, 2018, 4:15:39 PM6/2/18
to ipi-users
Hello,
I wrote a toy driver to see if I was able to interface a c++ driver with i-pi. There is only one file:
*************************************************
#include <cmath>
#include <cassert>
#include <iomanip>
#include <iostream>
#include <fstream>
#include <cstring>
#include <stdexcept>
#include "../external/sockets.h"
const int LENMSG = 12;
const int LENINIT = 2048;
const std::string STATUS = "STATUS ";
const std::string NEEDINIT = "NEEDINIT ";
const std::string HAVEDATA = "HAVEDATA ";
const std::string READY = "READY ";
const std::string GETFORCE = "GETFORCE ";
const std::string POSDATA = "POSDATA ";
const std::string FORCEREADY = "FORCEREADY ";
////////////////////////////////////////////////////////////////////////////////
int main(int argc, char** argv)
{
// Initialize defaults
int socket = 0;
int inet = 0;
int port = 31415;
std::string host = "localhost";
open_socket(&socket, &inet, &port, host.data());
// Variables needed for MD loop
char init_buffer[LENINIT + 1];
char header[LENMSG + 1];
std::vector<double> box = {100, 0, 0, 0, 100, 0, 0, 0, 100};
std::vector<double> boxi = {100, 0, 0, 0, 100, 0, 0, 0, 100};
std::vector<double> virial = {0, 0, 0, 0, 0, 0, 0, 0, 0};
int rid;
int buffl = LENMSG;
bool isinit = false;
bool hasdata = false;
int nat = 1;
double energy = 0.0;
int bsize = 0;
std::vector<double> buffer;
try {
while (true) {
while (true) {
buffl = LENMSG;
if (!hasdata) {
readbuffer(socket,header,buffl);
} else {
throw std::runtime_error("Wrapper did not ask for data yet");
}
if (strncmp(header,STATUS.data(),STATUS.length()) == 0) {
writebuffer(socket,READY.data(),READY.length());
} else {
break;
}
}
if (strncmp(header,POSDATA.data(),POSDATA.length()) == 0) {
readbuffer(socket, (char*) box.data(), 9*sizeof(double));
readbuffer(socket, (char*) boxi.data(), 9*sizeof(double));
readbuffer(socket, (char*) &nat, sizeof(int));
if (bsize == 0) {
bsize = 3*nat;
buffer = std::vector<double>(3*nat);
} else if (bsize != 3*nat) {
throw std::runtime_error("Number of atoms has changed.");
}
readbuffer(socket, (char*) buffer.data(), bsize*sizeof(double));
} else {
throw std::runtime_error("Wrapper did not send the positions");
}
std::cout << "x = " << buffer[0] << std::endl;
// Get forces here and store them in buffer
// Forces and energy and everything is input/output in au
// Test harmonic in X
double x = buffer[0];
energy = x*x;
buffer[0] = -2.0 * x;
buffer[1] = 0.0;
buffer[2] = 0.0;
virial[0] = buffer[0] * x;
//std::fill(buffer.begin(), buffer.end(), 0.0);
hasdata = true;
header[0] = '\0';
while (true) {
buffl = LENMSG;
if (hasdata) {
readbuffer(socket,header,buffl);
} else {
throw std::runtime_error("No data to sent to wrapper");
}
if (strncmp(header,STATUS.data(),STATUS.length()) == 0) {
writebuffer(socket,HAVEDATA.data(),HAVEDATA.length());
} else {
break;
}
}
if (strncmp(header,GETFORCE.data(),GETFORCE.length()) == 0) {
writebuffer(socket, FORCEREADY.data(), FORCEREADY.length());
writebuffer(socket,(char*) &energy,sizeof(double));
writebuffer(socket,(char*) &nat,sizeof(int));
writebuffer(socket,(char*) buffer.data(), bsize*sizeof(double));
writebuffer(socket,(char*) virial.data(),9*sizeof(double));
nat=0;
writebuffer(socket,(char*) &nat,sizeof(int));
} else {
throw std::runtime_error("Wrapper did not ask for forces");
}
hasdata = false;
}
} catch (const std::exception& e) {
std::cerr << " ** Error ** : " << e.what() << std::endl;
return 1;
}
return 0;
}
*************************
Then I take the init.pdb and input.xml from i-pi/examples/harmonic, and modify the address from harmonic to hostname:
**********************
<simulation verbosity='high'>
<output prefix='simulation'>
<properties stride='5' filename='out'> [ step, time{picosecond}, conserved{kelvin}, temperature{kelvin}, kinetic_cv{kelvin}, potential{kelvin}, pressure_cv{megapascal}] </properties>
<trajectory filename='pos' stride='5' format='xyz' cell_units='angstrom'> positions{angstrom} </trajectory>
<trajectory filename='force' stride='5' format='xyz' cell_units='angstrom'> forces{piconewton} </trajectory>
</output>
<total_steps>1000</total_steps>
<prng>
<seed>3348</seed>
</prng>
<ffsocket mode='unix' name='driver'>
<address>localhost</address>
</ffsocket>
<system>
<initialize nbeads='4'>
<file mode='pdb'> init.pdb </file>
<velocities mode='thermal' units='kelvin'>3683.412077</velocities>
</initialize>
<forces>
<force forcefield='driver'/>
</forces>
<ensemble>
<temperature units='kelvin'>3683.412077</temperature>
</ensemble>
<motion mode='dynamics'>
<dynamics mode='nvt'>
<thermostat mode='langevin'>
<tau units='femtosecond'>25</tau>
</thermostat>
<timestep units='femtosecond'>0.003</timestep>
</dynamics>
</motion>
</system>
</simulation>
****************************************************
This input works, and the simulation seems to update properly coordinates and forces. However, changing the number of beads to 1, makes i-pi not update the positions at all.
I am missing probably some information in the xml, but I haven't been able to find it.
Any hint/advice on what is going on will be greatly appreciated.
Thanks!
Marc
PS: the code has several unused variables and is not properly optimized. Sorry for that.
venkat kapil
Jun 2, 2018, 5:04:20 PM6/2/18
to ipi-users
The centre of mass is fixed by default. For a single particle it implies no movement. You can set it to false by adding <fixcom> False </fixcom> inside
<motion mode='dynamics'> ... </motion>
<motion mode='dynamics'> ... </motion>
Reply all
Reply to author
Forward
0 new messages