Write full box to topology file + script representation

Question

I would like topo writelammpsdata to print the full box, while it seems only to be printing one molecule.

I have loaded several molecules in a box with VMD - here working with water and benzene as an example. Although the box looks alright, I am not able to correctly topo writelammpsdata. Essentially, while I have N molecules, the lt file produced only contains 1 molecule.

I paste the whole script and the initial files. Please note that I am translating each molecule and picking alternatively a water or a benzene molecule. Also, is there a way to script a representation in VMD? Say I want each water to be licorice.

set side 50
set nmol_per_side 2

## Replicate the molecule

mol new water.pdb
set mol_idx 0

for {set i 0} {$i <= $nmol_per_side} {incr i} {
  for {set j 0} {$j <= $nmol_per_side} {incr j} {
    for {set k 0} {$k <= $nmol_per_side} {incr k} {

        set dx [expr $side/$nmol_per_side*$i]
    set dy [expr $side/$nmol_per_side*$j]
        set dz [expr $side/$nmol_per_side*$k]

        if { ($dx != 0) || ($dy != 0) || ($dz != 0) } {
        incr mol_idx
        puts "dx = $dx, dy = $dy, dz = $dz, mol_idx = $mol_idx"


            if { ($mol_idx % 2)==0 } {
           mol new water.pdb
        } else {
           mol new benzene.pdb
           }
        set molecule [atomselect $mol_idx "all"]
            $molecule moveby [ list $dx $dy $dz ]

            }
      }
    }
  }

# Create a cubic box
pbc set [ list $side $side $side ]
pbc box

topo retypebonds
topo guessangles
topo guessdihedrals

topo writelammpsdata full_box.lt full

You can really use any molecule. Here I am using water:

ATOM      1  O   HOH     1       0.000   0.000   0.000  1.00  0.00           O  
ATOM      2  H1  HOH     1       0.957   0.000   0.000  1.00  0.00           H  
ATOM      3  H2  HOH     1      -0.239   0.927   0.000  1.00  0.00           H  
CONECT    1    2    3  
CONECT    2    1  
CONECT    3    1  

and benzene:

ATOM      1  C   BEN     1       0.000   1.402   0.000  1.00  0.00           C  
ATOM      2  C   BEN     1       1.214   0.701   0.000  1.00  0.00           C  
ATOM      3  C   BEN     1       1.214  -0.701   0.000  1.00  0.00           C  
ATOM      4  C   BEN     1       0.000  -1.402   0.000  1.00  0.00           C  
ATOM      5  C   BEN     1      -1.214  -0.701   0.000  1.00  0.00           C  
ATOM      6  C   BEN     1      -1.214   0.701   0.000  1.00  0.00           C  
ATOM      7  H   BEN     1       0.000   2.494   0.000  1.00  0.00           H  
ATOM      8  H   BEN     1       2.156   1.247   0.000  1.00  0.00           H  
ATOM      9  H   BEN     1       2.156  -1.247   0.000  1.00  0.00           H  
ATOM     10  H   BEN     1       0.000  -2.788   0.000  1.00  0.00           H  
ATOM     11  H   BEN     1      -2.156  -1.247   0.000  1.00  0.00           H  
ATOM     12  H   BEN     1      -2.156   1.247   0.000  1.00  0.00           H  
CONECT    1    2    6  
CONECT    2    1    3  
CONECT    3    2    4  
CONECT    4    3    5  
CONECT    5    4    6  
CONECT    6    1    5  
CONECT    1    7  
CONECT    2    8  
CONECT    3    9  
CONECT    4   10  
CONECT    5   11  
CONECT    6   12  

Thank you

Marco

Answer 1

As per documentation (search for "molid"), TopoTools only works on one "VMD molecule" at a time, and by default this is just mol top. Your approach adds each new molecule to a new VMD molecule, and thus writelammpsdata writes only the last molecule you've created.

In your situation, the workflow I would adopt is:

1. Use writelammpsdata to convert water.pdb and benzene.pdb each to separate LAMMPS data files.

2. Read them both into a LAMMPS script using successive read_data commands (with the needed keywords, especially extra/* and offset).

3. In the same script, use replicate to multiply the water and benzene to the necessary quantity.

4. Equilibrate well (possibly at high temperature), adjusting the box size manually as desired.

Ordinarily I would recommend PACKMOL, which will combine water.pdb and benzene.pdb into (for example) 90water+70benzene.pdb, with some randomisation of coordinates to boot. You could then load that into VMD as a single molecule and topo writelammpsdata into a LAMMPS data file. But, for benzene, I imagine you need some care defining the dihedrals and impropers, and you might prefer to replicate in LAMMPS and know that all the benzenes have the same dihedral representation, rather than having TopoTools independently convert each one from PDB (or VMD's internal representation).