Dear Sayyed,

Speaking purely from C++ point of view, I think the application should take more care about array boundaries:

ifstream input1("in1.txt");

string atom[3];

double x[3], y[3], z[3];

int i = 0;

while (input1 >> atom[i]) //read solvent input file

{

input1 >> x[i];

input1 >> y[i];

input1 >> z[i];

i++;

}

like here - what would happen if we have a large in1.txt with more elements than array can hold? And what if the format is not correct - that would lead to errors, that are not checked for.

And why would you need multiple local arrays of size 30000 - I think you can utilize memory in a much better way if you use sets or vectors with the right dimension.

Best regards,

 Stan

Good tips stan, I'm beginner with c++, that's why I asked for your comments.

Cheers

@Sayyed,

There is nothing wrong with writing your own code for wanting to solve these kind of problems. However, you should notice that there already exist codes for solving the same problems in a more optimized way. Randomly positioning the solvent molecules around solute is not the most optimized way of doing it, since a long protocol will be needed to equilibrate solvent around solute. Another suggestion is related to the existence of other species in the solvent, for example, ions.

Some examples of existing codes are listed in the following:

1) 'solvate' module of VMD program, which attempts to place every water molecule in a position for which the potential energy of system has a minimum value, by using Monte Carlo simulation.

2) charmm program, where you can run scripts to use already equilibrated water box (which is easy to find out), and then try to insert the solute into it, by removing these water molecules which overlap with solute atoms.

Dear  Hiqmat, I already know that there are more sufficient programs and tools for solvation. I just did it as a training as I said I'm beginner with c++. Besides those tools you said, lammps has a great module using grand canonical monte carlo for solvaton.

Cheers.

Hi Sayyed,

I would recommend, as start, to also include a few Monte Carlo steps of rotations of solvent molecule after choosing a random position (as already implemented in your code). These steps could help to find the best orientation of the solvent molecule.

I'm working on it, thank you so much. I also thinking about including a simple non bonded inter-atomic potential for improvement.

Sayyed,

Please have a look at Packmol. http://www.ime.unicamp.br/~martinez/packmol/home.shtml

It is a wonderful program. I have used this program extensively for my computer simulations. 

Thank you Jestin.

You may also take a look at  conventional optimized solvation plugins  like  the one with VMD to  improve it.