How to prepare a solvated system of large molecules automatically

The Packmol distribution includes the solvate.tcl script, which is used to solvate large molecules, usually proteins, with water and ions (Na+ and Cl-).  Given the PDB file of the biomolecule, just run the script with:
        solvate.tcl PROTEIN.pdb

And the script will create a input file for packmol called packmol_input.inp. With this file, run Packmol with:
        packmol < packmol_input.inp

And your large molecule will be solvated by a shell of 15. Angs. of water, and ions to keep the system neutral and a physiological NaCl concentration of 0.16M. The script usually makes reasonable choices for every parameter (number of water molecules, number of ions, etc.), but these may be controlled manually with additional options, as described below:
        solvate.tcl structure.pdb -noions
        solvate.tcl structure.pdb -shell 15.  -charge +5  -density 1.0 -i pack.inp -o solvated.pdb

Where: structure.pdb is the pdb file to be solvated (usually a protein)

"15." is the size of the solvation shell. This is an optional parameter. If not set, 15. will be used.

"-charge +5" is the total charge of the system, to be neutralized. This is also and optional parameter, if not used, the package considers histidine residues as neutral, Arg and Lys as +1 and Glu and Asp as -1. The Na+ and Cl- concentrations are set the closest possible to 0.16M, approximately the physiological concentration.  Alternatively, use the -noions to not add any ions, just water.

1.0 is the desired density. Optional. If not set, the density will be set to 1.0 g/ml.

solvated.pdb: is the (optional) name for the solvated system output file. If this argument is not provided, it will be the default solvated.pdb file.

pack.inp: is the (optional) name for the packmol input file that will be generated. If not provided, packmol_input.inp will be used.

All these options are output when running the "solvate.tcl" script without any parameter. The script also outputs the size of the box and the suggested periodic boundary condition dimensions to be used.


obtained from Packmol user-guide (http://www.ime.unicamp.br/~martinez/packmol/userguide.shtml)

How to compile in Cygwin console

To compile in Cygwin console, additional installation is often required.
For example, gfortran is required to compile (build) Packmol.
However, default cygwin installation process does not include compilers (e.g., gcc, gfortran).

When installing Cygwin, we need to select necessary packages.
In a package search window, type "gcc" and then gcc packages will appear.
Change the default option into "install" or we can select particular packages in a subfolder.

In this way, other packages can be installed.
For example, to prepare MD simulation environment on Windows OS system, we will need:
required package of cygwin + python (for Moltemplate) + gcc (for Packmol) + make (for Packmol)
(gcc package includes gfortran)

How to use Moltemplate on Windows OS systems

Moltemplate (http://www.moltemplate.org/) is a powerful tool for configuring molecules and assigning force field parameters.
However, its use on Windows OS is somewhat tricky.

1. Install Cygwin (https://cygwin.com/install.html).  Installation of required packages + python should be good enough.

2. Copy all Moltemplate files to C:\cygwin64\home\moltemplate

3. Open "C:\cygwin64\home\username\.bash_profile" in a text editor.  Add a new PATH (in yellow) to a section "# Set PATH so it includes user's private bin if it exists".  For example,

# Set PATH so it includes user's private bin if it exists
# if [ -d "${HOME}/bin" ] ; then
   PATH="/bin/:/home/moltemplate/"
# fi

In Cygwin console (terminal), we can check if the PATH is correctly set by typing "echo $PATH"

5. Place examples (http://www.moltemplate.org/) in C:\cygwin64\home\moltemplate\examples

6. In Cygwin command window (console or teminal), navigate an example folder
(e.g., cd /home/moltemplate/examples/all_atom_examples/force_field_explicit/waterSPCE+Na+Cl/moltemplate_files)

7. In Cygwin command window, type:

moltemplate.sh  -atomstyle "full"  system.lt
or
moltemplate.sh  -xyz  system.xyz  -atomstyle "full"  system.lt
(This will use a premade xyz configuration of system, e.g., system.xyz).

LAMMPS input and data files will be created.  (e.g., system.in and system.data)

LAMMPS - Parallel computing on a multi-core Windows flatform

LAMMPS, a molecular dynamics simulator supports parallel computing.
To use such parallel function, we need to install MPI. (MPICH2, openMPI, or MS-MPI)


For Windows 7,
MPICH2 can be downloaded at  http://www.mpich.org/downloads/
openMPI can be downloaded at  https://www.open-mpi.org/software/
MS-MPI can be downloaded at  http://www.microsoft.com/en-au/download/details.aspx?id=36045


After installing,
1. two excutable files (mpiexec.exe and smpd.exe) need to be copied to LAMMPS bin folder (with lmp_serial.exe and lmp_mpi.exe).
2. open a command prompt (cmd.exe) as an administrator and then type the following commands in sequence.
     smpd -install
     mpiexec -remove
     mpiexec -register     (set up "username" and "password")
     mpiexec -validate     (it should return "SUCCESS")
     smpd -status     (it should return "smpd running on 'hostname')


Now, we are ready to run parallel computing.
When running, "-localonly '# of processors' " option should be used.
mpiexec -localonly 4 lmp_mpi < example.in

To make sure, just compare results of both serial and parallel simulations using same input script.
e.g.,
lmp_serial < example.in
mpiexec -localonly 4 lmp_mpi < example.in

[C/C++] read xyz file and convert element symbol into atomic number

// This program reads a text file containing molecular configuration (element symbol, and xyz coordinates), and convert element symbol into atomic number.
 
#include <stdio.h>
 
#define max 1000
 
    int element[max], t_element=0, c_element[max];
    float x[max], y[max], z[max], t_x, t_y, t_z, c_x[max], c_y[max], c_z[max];
    int total_n=0;
    int count_atom=0;
    void list_atoms(int count_atom, int t_element, float *t_x, float *t_y, float *t_z);
 
 
void main(void)
{
FILE *xyzfile = fopen("example.txt", "r");
    if(xyzfile==NULL) {
        printf("\nError: Unable to open file.\n");
        exit(1);
    } else {
        fscanf(xyzfile, "%d\n", &total_n);
//       printf("\n\tTotal number of atoms : %d\n", total_n);
        char description[100];
        fgets(description, 100, xyzfile);
//       printf("\nDescription : %s\n", description);
 
        char e_type[2];
        int i=1;           // [1] is for first atom (third line)
        while(fscanf(xyzfile, "%s %f %f %f\n", e_type, &x[i], &y[i], &z[i]) != EOF){       // loop through and store elemental symbol, x, y, and z values into the array
 
            if (strcmp(e_type, "H")==0 || strcmp(e_type, "h")==0) {                        // convert elemental symbol into elemetal number
                element[i] = 1;
                t_element=element[i];
                t_x=x[i]; t_y=y[i]; t_z=z[i];
                count_atom++;
                list_atoms(count_atom, t_element, &t_x, &t_y, &t_z);
                i++;
            } else if(strcmp(e_type, "C")==0 || strcmp(e_type, "c")==0) {
                element[i] = 6;
                t_element=element[i];
                t_x=x[i]; t_y=y[i]; t_z=z[i];
                count_atom++;
                list_atoms(count_atom, t_element, &t_x, &t_y, &t_z);
                i++;
            } else if(strcmp(e_type, "N")==0 || strcmp(e_type, "n")==0) {
                element[i] = 7;
                t_element=element[i];
                t_x=x[i]; t_y=y[i]; t_z=z[i];
                count_atom++;
                list_atoms(count_atom, t_element, &t_x, &t_y, &t_z);
                i++;
            } else if(strcmp(e_type, "O")==0 || strcmp(e_type, "o")==0) {
                element[i] = 8;
                t_element=element[i];
                t_x=x[i]; t_y=y[i]; t_z=z[i];
                count_atom++;
                list_atoms(count_atom, t_element, &t_x, &t_y, &t_z);
                i++;
            } else if(strcmp(e_type, "Na")==0 || strcmp(e_type, "na")==0 || strcmp(e_type, "NA")==0) {
                element[i] = 11;
                t_element=element[i];
                t_x=x[i]; t_y=y[i]; t_z=z[i];
                count_atom++;
                list_atoms(count_atom, t_element, &t_x, &t_y, &t_z);
                i++;
            } else if(strcmp(e_type, "Mg")==0 || strcmp(e_type, "mg")==0 || strcmp(e_type, "MG")==0) {
                element[i] = 12;
                t_element=element[i];
                t_x=x[i]; t_y=y[i]; t_z=z[i];
                count_atom++;
                list_atoms(count_atom, t_element, &t_x, &t_y, &t_z);
                i++;
            } else if(strcmp(e_type, "Ca")==0 || strcmp(e_type, "ca")==0 || strcmp(e_type, "CA")==0) {
                element[i] = 20;
                t_element=element[i];
                t_x=x[i]; t_y=y[i]; t_z=z[i];
                count_atom++;
                list_atoms(count_atom, t_element, &t_x, &t_y, &t_z);
                i++;
            } else {
                element[i] = 0;
            }
    }
    fclose(xyzfile);
 
    printf("\nNumber of atoms : %d", count_atom);
}
}
 
 
void list_atoms(int count_atom, int t_element, float *t_x, float *t_y, float *t_z)
{
    int k = count_atom;
 
    c_element[k] = t_element;
    c_x[k] = *t_x;
    c_y[k] = *t_y;
    c_z[k] = *t_z;
    printf("\nafter function, converted element %-3d = %-2d %.5f %.5f %.5f", count_atom, c_element[k], c_x[k], c_y[k], c_z[k]);
}
 

An easy way to create initial molecular system for computer simulation

Creating a single molecular structure is relatively easy if we use a molecule editor (e.g, Avogadro, Materials Studio, etc).

However, a molecular system consisting of many molecules cannot be easily prepared by those editors.

"Packmol" is a good software for readily creating such a large system.

http://www.ime.unicamp.br/~martinez/packmol/

or  http://leandro.iqm.unicamp.br/packmol/versionhistory/

1. Install "Packmol" as instructed. (it can easily be built in Cygwin64 console.)

2. Prepare component molecules in pdb, xyz, etc file format.

    (component molecules can be DPD structure.)

3. Make an input script 

4. Run "packmol < filename.inp"

An example Packmol input script is below:

# Every atom from diferent molecules will be far from each other at least 2.0 Angstroms in the system.

tolerance 2.0 

# Coordinate file types will be in pdb format (keyword not required for pdb file format, but required for tinker, xyz or moldy).

filetype xyz

# The output configuration filename

output 100water_5Na_5Cl.xyz

# The first three numbers are the minimum x, y, z coordinates for this molecules, the last three are maximum coordinates. 

structure Na.xyz

  number 5

  inside box 0. 0. 0. 20. 20. 20.

end structure

structure Cl.xyz

  number 5

  inside box 0. 0. 0. 20. 20. 20.

end structure 

structure water.xyz 

  number 100

  inside box 0. 0. 0. 20. 20. 20.
# inside box  Xmin Ymin Zmin  Xmax Ymax Zmax

# inside sphere  10.0 10.0 15.0 8.0
# inside sphere  a b c d  (a,b,c,d in a sphere equation, (x-a)^2 + (y-b)^2 + (z-c)^2 - d^2 = 0,  d is the diameter)
# inside cylinder  0.0 0.0 0.0  1.0 0.0 0.0  5.0  10.0    
# inside cylinder  x1 y1 z1  x2 y2 z2  d  l  (direction_of_line_in_vector_(x1 y1 z1  x2 y2 z2),  distance_to_the_line_(d),  length_of_cylinder_(l))

end structure