Doubts with bonds_type

[ProfLeao]

Hi Andrew. I've been looking at scripts (sh) but I'm not very good with this programming language, I didn't understand it very well. I'm in doubt about the bonds_type property of the GPSettings class, in genpoly. I know that genpoly_lt does not check the bond and atom type, but moltemplate.sh does. I'm trying to use types present in force field, like "084_084" from oplsaa, but moltemplate says the type doesn't have the necessary parameters. Where can I learn more about the valid types of bonds_type and bonds_atoms?

[jewettaij]

Hi Andrew. I've been looking at scripts (sh) but I'm not very good with this programming language, I didn't understand it very well. I'm in doubt about the bonds_type property of the GPSettings class, in genpoly. I know that genpoly_lt does not check the bond and atom type, but moltemplate.sh does.

Hi Reginaldo I think I understand what might be causing your confusion. In all of the moltemplate examples which use OPLSAA, the user is not required to specify the type of bond connecting pairs of atoms. The user only has to specify which atoms are bonded together (placing them in the write("Data Bonds List") section). Then moltemplate will lookup the appropriate @bond type using the rules in the "Data Bonds By Type" section of the "oplsaa.lt" file. Unfortunately "genpoly_lt.py" is currently not smart enough to do this. When using "genpoly_lt.py", the user must explicitly specify the @bond type for bonds between atoms in different monomers (such as "084_084").

It's important that you bring this to my attention, because I had forgotten this issue. I have only used "genpoly_lt.py" with coarse-grained polymers (which don't use OPLSAA), so I did not think about how confusing this must be.

But it sounds like you are attempting to do the correct thing: You are manually specifying the bond types ("084_084"). That's good. I don't know why moltemplate.sh is complaining...

I'm trying to use types present in force field, like "084_084" from oplsaa, but moltemplate says the type doesn't have the necessary parameters. Where can I learn more about the valid types of bonds_type and bonds_atoms?

From your description, it sounds like the problem becomes apparent when you run "moltemplate.sh" on your LT files. (It sounds like "genpoly_lt.py" runs without error messages.) Is this correct? If you haven't already, can you send me your LT files? (Email is okay.) Please also explain how you use genpoly_lt.py to create your LT files.

[ProfLeao]

Hi Andrew, thanks for the feedback. I made a script to automate what I've been doing here, you can easily reproduce it and then check for the error. I sent you by email.

 Syntax error: Missing bond coeff.

  No coeffs for the "@/bond:_sub1/084_084" bond type have been
defined, but a reference to that bond type was discovered
near "HDPE_GP8553_CADS8.LT", line 42799.   Check this file and also check
your "bond_coeff" commands or your "Data Bond Coeffs" section.
---------------------------------------------------------------------
(To continue anyway, run moltemplate using the "-nocheck" argument.)

[ProfLeao]

To leave it documented for other users I will be more specific. I'm using the following groups:

import "oplsaa.lt"
import "loplsaa.lt" 

CH2 inherits OPLSAA {

  write("Data Atoms") {
    $atom:C  $mol:... @atom:81L    0.0  0.000000     0.000000      0.000000
    $atom:H1 $mol:... @atom:85LCH2 0.0  0.000000     0.631044      0.892431
    $atom:H2 $mol:... @atom:85LCH2 0.0  0.000000     0.631044     -0.892431
  }

  write('Data Bond List') {
    $bond:CH1 $atom:C $atom:H1
    $bond:CH2 $atom:C $atom:H2
  }

import "oplsaa.lt"
import "loplsaa.lt" 

CH3 inherits OPLSAA {

  write("Data Atoms") {
    $atom:C  $mol:... @atom:80L    0.0  0.000000     0.000000      0.000000
    $atom:H1 $mol:... @atom:85LCH3 0.0  0.000000     0.631044      0.892431
    $atom:H2 $mol:... @atom:85LCH3 0.0  0.000000     0.631044     -0.892431
    $atom:H3 $mol:... @atom:85LCH3 0.0 -0.892431    -0.631044      0.000000
  }

  write('Data Bond List') {
    $bond:CH1 $atom:C $atom:H1
    $bond:CH2 $atom:C $atom:H2
    $bond:CH3 $atom:C $atom:H3
  }

A function to automate the generation of cuts and sequences.

def sequence(mon_lst, pol_deg, dims_lattice, modo="file", random= False):
    try:
        num_sitios = dims_lattice[0] * dims_lattice[1] * dims_lattice[2]
        num_cadeias = num_sitios // pol_deg
        pol_cap, pol_body = mon_lst[0], mon_lst[1]
    except:
        return ValueError

    with open("cuts.txt", 'w') as cuts: 
        line = 0
        with open("sequence.txt", 'w') as seq:
            for i in range(num_cadeias):
                seq.write(pol_cap + "\n")
                line+=1
                for j in range(pol_deg-2):
                    seq.write(pol_body + "\n")
                    line+=1
                seq.write(pol_cap + ".rot(180,0,0,1)" + "\n")
                line+=1
                if i != (num_cadeias - 1):
                    cuts.write(f"{line}\n")
                else:
                    pass

def main():
    mon_lst = sys.argv[1:3]
    sys.argv[3:7] = [int(i) for i in sys.argv[3:7]]
    pol_deg, dims_lattice =  sys.argv[3], sys.argv[4:7]
    try:
        sequence(mon_lst, pol_deg, dims_lattice)
    except ValueError:
        print("**** ERRO ****\n", "Parâmetros inválidos!")
        exit()
    print("Sequencias de monômeros e de cortes geradas com sucesso!")
if __name__ == "__main__":
    main()

And a routine for automating the polyethylene generation procedure.

from PolymerCpp.helpers import getCppWLC
import sequence_generator as sqg
import numpy as np
import sys
from new_genpoly_lt import GenPoly
import subprocess
import os

def main():
    """
        O programa recebe o grau de polimerização e o número de cadeias a serem 
        simuladas, diretamente por linha de comando. 

        The script receives the degree of polymerization and the number of chains to be
        simulated, directly from the command line.
    """

    polim_degree, num_chains = int(sys.argv[1]),int( sys.argv[2])
    pers_length = 2
    ver_pers_len = input(
        f"Using persistence length equal to  {pers_length}. " +
        "Do you want to modify this parameter? [y/n]: "
    )
    if ver_pers_len in ['y', 'Y']:
        pers_length = int(
            input(
                "Enter a new persistence length" +
                "int: "
            )
        )
    chains = getCppWLC(polim_degree *num_chains, pers_length)

    sqg.sequence(
        ["CH3","CH2"], 
        polim_degree,
        [polim_degree *num_chains,1,1]
    )

    max_x, min_x = np.amax(chains[:,0]), np.amin(chains[:,0])
    max_y, min_y = np.amax(chains[:,1]), np.amin(chains[:,1])
    max_z, min_z = np.amax(chains[:,2]), np.amin(chains[:,2])
    gp = GenPoly()
    gp.coords_multi = list(np.split(chains[:-1],num_chains))

    gp.ParseArgs(
        [
            '-polymer-name', 'HDPE2911',
            '-monomer-name', 'alcano',
            '-sequence', 'sequence.txt',
        '-header', 'import "oplsaa.lt"',
            '-header', 'import "loplsaa.lt"',
            '-header', 'import "ch2group.lt"',
        '-header', 'import "ch3group.lt"',
        '-cuts', 'cuts.txt', 
            '-bond', '??', '??', '??',
        ]
    )
    with open(
        f"HDPE_GP{polim_degree}_CADS{num_chains}.LT", 
        'w'
    ) as molt_file:
        gp.WriteLTFile(molt_file)
        molt_file.write(
            'write_once("Data Boundary") {\n' +\
            f"    {1.1 * min_x}  {1.1 * max_x}  xlo xhi\n" +\
            f"    {1.1 * min_y}  {1.1 * max_y}  ylo yhi\n" +\
            f"    {1.1 * min_z}  {1.1 * max_z}  zlo zhi\n"
            "}"
        )
    return True
if __name__ == "__main__":
    main()
    print("Input successfully generated!")

PolymerCPP library can be installed with pip install PolymerCpp.

The line chains = getCppWLC(polim_degree *num_chains, pers_length) generates an 1A spaced chain. I think it is not necessary to take care of that since I use minimize in Lammps. Right?

Thanks.

[jewettaij]

I apologize for taking so much time to reply. I should have warned you that I am very busy. But I have your files and I will work on it by the week end.

[jewettaij]

Hi Reginaldo My apologies for taking so long to reply. My career has changed and I cannot spend as much time on moltemplate as I did in the past. But I would like to solve your problem so that other people do not struggle as you did. Unfortunately, I cannot get your code to run. Perhaps I am running it the wrong way. When I run your code this way:

./my_auto_pol_gen.py

...I get the following error:

Traceback (most recent call last):
  File "./my_auto_pol_gen.py", line 3, in <module>
    from PolymerCpp.helpers import getCppWLC
  File "/home/jewett/Downloads/moltemplate_question_2021-11-06__reginaldo.junior@ifmg.edu.br/venv_PolymerCpp/lib/python3.8/site-packages/PolymerCpp/helpers.py", line 9, in <module>
    import PolymerCppCore
ImportError: /home/jewett/Downloads/moltemplate_question_2021-11-06__reginaldo.junior@ifmg.edu.br/venv_PolymerCpp/lib/python3.8/site-packages/PolymerCppCore.cpython-38-x86_64-linux-gnu.so: undefined symbol: _ZSt28__throw_bad_array_new_lengthv

Am I using your code the correct way?

If it matters, I installed the dependencies using:

pip install PolymerCpp
pip install moltemplate

If too much time has gone by and you are no longer interested in this issue, I will not be offended.

[bond]

Hi, @bond tags me in the thread.

On Sun, Nov 7, 2021, 00:26 Andrew Jewett @.***> wrote:

Hi Andrew. I've been looking at scripts (sh) but I'm not very good with this programming language, I didn't understand it very well. I'm in doubt about the bonds_type property of the GPSettings class, in genpoly. I know that genpoly_lt does not check the bond and atom type, but moltemplate.sh does.

Hi Reginaldo I think I understand what might be causing your confusion. In all of the moltemplate examples which use OPLSAA, the user is not required to specify the type of bond connecting pairs of atoms. The user only has to specify which atoms are bonded together (placing them in the write("Data Bonds List") section). Then moltemplate will lookup the appropriate @bond https://github.com/bond type using the rules in the "Data Bonds By Type" section of the "oplsaa.lt" file. Unfortunately "genpoly_lt.py" is currently not smart enough to do this. When using "genpoly_lt.py", the user must explicitly specify the @bond https://github.com/bond type for bonds between atoms in different monomers (such as "084_084").

It's important that you bring this to my attention, because I had forgotten this issue. I have only used "genpoly_lt.py" with coarse-grained polymers (which don't use OPLSAA), so I did not think about how confusing this must be.

But it sounds like you are attempting to do the correct thing: You are manually specifying the bond types ("084_084"). That's good. I don't know why moltemplate.sh is complaining...

I'm trying to use types present in force field, like "084_084" from oplsaa, but moltemplate says the type doesn't have the necessary parameters. Where can I learn more about the valid types of bonds_type and bonds_atoms?

From your description, it sounds like the problem becomes apparent when you run "moltemplate.sh" on your LT files. (It sounds like "genpoly_lt.py" runs without error messages.) Is this correct? If you haven't already, can you send me your LT files (email is okay)? Please also explain how you use genpoly_lt.py to create your LT files.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/jewettaij/moltemplate/issues/62#issuecomment-962523318, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAABNMYOYQT4YTBNM6K4DPDUKXBSNANCNFSM5HO6PAQA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[jewettaij]

Hi, @bond tags me in the thread. […] My original response: "Lucky you. You can learn all about how to simulate polymers in LAMMPS by reading this thread. Cheers. -A"

I did not mean to sound so flippant. If you know how to prevent receiving these messages, let me know. Either way, hopefully, we will get this issue resolved soon. -A 2021-12-28

[jewettaij]

Hi Reginaldo. Perhaps I can solve your problem. I just posted commit e37b1f5, which allows the user to omit the bond's type. Here is the documentation for this feature:

    -bond a1 a2
             Add a bond between successive monomers between atoms named a1, a2.
         Here we omit the bond's type.  Omitting the bond's type is allowed
             if you are using a force field (like OPLSAA), which looks up bond
             types according to rules defined in the force field.

You can also use the original syntax, where the bond type (eg "084_084") is explicitly specified by the user.

     -bond btype a1 a2

(In your case, btype="084_084".) I realize you have probably moved on to solve your problem in some other way, and I apologize for that. But I appreciate you posting your question because I suspect other people who use "genpoly_lt.py" will run into this problem too (since most of them also use force fields too). This will hopefully help them. Thanks. -A P.S. If I don't hear from you in a a couple weeks, I will close this issue.

[ProfLeao]

Hello, @jewettaij. I'm coming back from my vocations today, I'll check it tomorrow. Ok? Thanks for returning!

[jewettaij]

Hi Reginaldo. I think I'm going to close this issue. Feel free to open it again later if you have time to revisit the problem. I hope things are going well.

[ProfLeao]

Hi @jewettaij , thank you very much for the effort to help with this issue. I've had several health problems in the family in the last few months, so I've only now taken up the issue again. I'm testing everything now.

[jewettaij]

I am sorry to hear that. I've been living with some sick family members too. This kind of problem can consume you. Take care. I'll be here to look at this issue whenever you are ready.

[ProfLeao]

I am sorry to hear that. I've been living with some sick family members too. This kind of problem can consume you. Take care. I'll be here to look at this issue whenever you are ready.

Perfect! Genpoly is working exactly as expected. Thank you so much.

[jewettaij]

I am sorry to hear that. I've been living with some sick family members too. This kind of problem can consume you. Take care. I'll be here to look at this issue whenever you are ready.

Perfect! Genpoly is working exactly as expected. Thank you so much.

I'm relieved. Thanks for getting back to me, Reginaldo. Good health for you and your family!

Andrew