Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

LAMMPS DATA WRITER #3197

Closed
ksy141 opened this issue Mar 31, 2021 · 3 comments
Closed

LAMMPS DATA WRITER #3197

ksy141 opened this issue Mar 31, 2021 · 3 comments

Comments

@ksy141
Copy link

ksy141 commented Mar 31, 2021

Hello, I'm trying to write a LAMMPS data file from a coordinate file. Consider a hydrogen gas,

cat > h2.gro << EOF
MDAnalysis BEST!!
 2
    1HYDR   HH1    1   3.000   3.000   3.000
    1HYDR   HH2    2   3.000   3.000   3.500
    6.00000   6.00000   6.00000
EOF

Now let's make a data file with MDAnalysis. In order to make a LAMMPS data file, the atom types should be changed to positive integers.

import MDAnalysis as mda
u = mda.Universe('h2.gro')
u.atoms.types = ['1', '1']
u.atoms.write('h2.data')

However, I received the error from the LAMMPS.py (~line 431).

attrs = [('bond', 'bonds'), ('angle', 'angles'), ('dihedral', 'dihedrals'), ('improper', 'impropers')]
for btype, attr_name in attrs:
    features[btype] = atoms.__getattribute__(attr_name)
    self.f.write('{:>12d}  {}\n'.format(len(features[btype]), attr_name))
    features[btype] = features[btype].atomgroup_intersection(atoms, strict=True)

I think try, except should be added like this?

for btype, attr_name in attrs:
    try:
        features[btype] = atoms.__getattribute__(attr_name)
        self.f.write('{:>12d}  {}\n'.format(len(features[btype]), attr_name))
        features[btype] = features[btype].atomgroup_intersection(atoms, strict=True)
    except:
        pass

After adding try, except the above example should work. Now let's add a bond between the first hydrogen atom (index: 0) and second hydrogen atom (index: 1)!

import MDAnalysis as mda
u = mda.Universe('h2.gro')
u.atoms.types = ['1', '1']
u.add_TopologyAttr('bonds', [[0, 1]])
u.atoms.write('h2.data')

I received the following error:
TypeError: LAMMPS DATAWriter: Trying to write bond, but bond type ('1', '1') is not numerical.

Trying to change the bond type like this is not allowed, with the error message saying AttributeError: can't set attribute

u.bonds[0].type = 1

Thanks! Anyone can help me how to set a bond properly so that MDAnalysis can understand?

MDAnalysis 1.0.1

@ksy141
Copy link
Author

ksy141 commented Mar 31, 2021

How about the following changes in LAMMPS.py (~ line 340)? I confirmed this works for bonds and angles for a simple water system. Commented lines are the lines included in the original code.

def _write_bonds(self, bonds):
    self.f.write('\n')
    self.f.write('{}\n'.format(btype_sections[bonds.btype]))
    self.f.write('\n')

    bt = {}; t=0
    for bond, i in zip(bonds, range(1, len(bonds)+1)):
        if bond.type not in bt.keys():
            t += 1
            bt[bond.type] = t
        self.f.write('{:d} {:d} '.format(i, bt[bond.type])+\
                ' '.join((bond.atoms.indices + 1).astype(str))+'\n')
    
#        try:
#            self.f.write('{:d} {:d} '.format(i, int(bond.type))+\
#                    ' '.join((bond.atoms.indices + 1).astype(str))+'\n')
#        except TypeError:
#            raise_from(TypeError('LAMMPS DATAWriter: Trying to write bond, '
#                            'but bond type {} is not '
#                            'numerical.'.format(bond.type)),
#                        None)

@1ut
Copy link
Contributor

1ut commented Apr 2, 2021

The codes you provided have been ​rewritten in develop branch like this commit. (0d951a2#diff-b36aa1d8c19a92c6f8918f35d8c0ea354fd051e64b40e789452551e7bc0792ba , PR #2759)

Also, I think some issues you provided have been solved in develop branch.

You may try a development version. The method to build it is here. (https://userguide.mdanalysis.org/stable/contributing_code.html#build-mdanalysis-develop)

@ksy141
Copy link
Author

ksy141 commented Apr 2, 2021

Thanks, closing the issue.

@ksy141 ksy141 closed this as completed Apr 2, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants