choderalab / openmoltools

An open set of tools for automating tasks relating to small molecules
MIT License
63 stars 30 forks source link

packmol bond order error #287

Open guoweiqi opened 5 years ago

guoweiqi commented 5 years ago

Trying to solvate a small molecule from trajectories or straight from PDB files for whatever reason inputs bond orders that are greater than 3... Any idea why this might be happening? Here's a short block of code that is triggering the error:

n_monomers = [1,200]
mol_traj = md.load('guest-4/guest-4.mol2')
water_traj = md.load('guest-4/water.mol2')

mol_traj.save_pdb('guest-4/guest-4_vac.pdb')

box_size = packmol.approximate_volume_by_density(['OC1CCCC1', 'O'], n_monomers)

packed_traj = packmol.pack_box([mol_traj, water_traj], n_molecules_list=[1, 200], box_size=box_size)
packed_traj.save_pdb('guest-4/guest-4_solv.pdb')

Test files can be found here:

https://github.com/MobleyLab/benchmarksets/tree/master/input_files/cd-set1

Eli-VDB commented 4 years ago

I have run into a similar problem, and I solved it in the following way, perhaps good to incorporate something similar in packmol.py: I have manually adjusted the packmol.py file (in ~/anaconda3/envs/openmm/lib/python3.7/site-packages/openmoltools) as the bonds are defined as an array of 1 x 4; at line 196 insert the code snippet below. When loading the molecules the bond type and order should be in the 3rd and fourth column, this is however not the case (I have not investigated this into detail but just wrote a quick workaround). The solve this is, put to order and type to zero which corresponds to a bond order and type of None, Antechambers will normally put this correct based on the bonds.

I believe this can be done in an alternative way without having the chance of losing information in this way.

#cast atom type and order to 0;
   if len(bonds[0])==4:
        for idx,frame in enumerate(bonds):
              bonds[idx][2] = 0;
              bonds[idx][3] = 0;