Fix bond/create/angle form multiple bonds instead of one

LAMMPS LAMMPS General Discussion
1

Hello everyone!

I am trying to simulate the following reaction

BB-1-2 + 3-4 —> BB-1-2-3-4

where a bond is formed between 2-3. In the meantime, I wish to create a new angle between 1-2-3.
I tried fix bond/create/angle, and it gives 2 angles of the same type for the above reaction (1-2-3 & 2-3-4).
Since there are multiple 1-2 side chains in the backbone, I assume it would be really complex to use fix bond/react.
Is there any simple way I can do this? Any suggestion would be appreciated!

Thanks,
Yang

2

This is the expected behavior. When you create a bond connecting a chain, there will be two angles added. This is what is consistent with all molecular force fields that I know. What you are asking for is rather unusual and inconsistent.

This is a question for @jrgissing. I don’t know the fix well enough.

3

Hi!

Thank you very much for your reply!

I understand that 2 angles will form when 2 chains connect. I think my question can be easily solved if fix bond/create/angle can form different types of angles. However, this seems to fall to fix bond/react.

Anyway, thank you very much.

Best,
Yang

4

This fix was created for use with spring-bead polymers.

Creating complex bonded interactions is not a problem the fix tries to solve, and it is hellishly complicated to do so in parallel.

Yes, fix bond/react was devised to address such issues and looking at the source code for the fix confirms my assessment of the complexity of the task.

5

fix bond/react was designed to handle complex chemistries, including highly branched or cyclic structures, so the existence of side chains is not an issue

6

Hi Jacob,

After reading the fix bond/react I am still bit confused about the file format.

So I wish to simulate the following cross-linking reaction

BB-1-2 + 3 ----> BB-1-2-3
BB-1-2-3 + 2-1-BB ----> BB-1-2-3-2-1-BB

And I want an angle in for atoms 2-3-2.

I have a map for the step1

          3 equivalence
          2 edgeIDs
          
          InitiatorIDs
          
          2
          3
          
          EdgeIDs
          
          1
          3
          
          Equivalences
          
          1  1
          2  2
          3  3

And a file for unreacted molecule

              3 atoms
              1 bonds
      
      Coords
      
      1    0.00000   0.00000   0.00000
      2    0.30000   0.00000   0.00000
      3    0.83000   0.00000   0.00000
      
      Types
      
      1        3   
      2        4   
      3        5   
      
      Bonds
      
      1   2      1      2

File for reacted molecule

              3 atoms
              2 bonds

      Coords
      
      1    0.00000   0.00000   0.00000
      2    0.30000   0.00000   0.00000
      3    0.83000   0.00000   0.00000
      
      Types
      
      1        3   
      2        4   
      3        5   
      
      Bonds
      
      1   2      1      2
      2   3      2      3

The simulation is stuck when I invoke fix bond/react. I know I get something wrong but I just can’t spot. Could you help with it?

Thank you!

7

Can you provide the full input deck so I can take a closer look?

8

The page says new users cannot upload attachments. I have pasted the input below. The data file is too large to be pasted here. Do you also need it?

Type 1,2 are BB,
Type 3,4 are sidechains
Type 5,6 are in a rigid body and type 6 do not interact.

I want type 4 and type 5 to form bonds.

units           lj
dimension       3
boundary        p p p

atom_style      full 

read_data       one_fiber.data extra/bond/per/atom 32 extra/special/per/atom 32 extra/angle/per/atom 32 extra/dihedral/per/atom 32 

variable        tEqui equal  4000
variable        tRun  equal  25000000
variable        tprod equal  400000
variable        deltat equal  0.01
#=========== frequency for bond creation =================
variable        Nevery equal  400

#=========== Parameters of the system ====================
variable        temp   equal 1.0
variable        T_mix  equal  ${temp}
#============ dielectric of the system =====================
variable        epsilon_r equal 1
#================ diameter size particle ====================
variable        sigma_a  equal  0.1  # filler/ligand/linker beads 
variable        sigma_b  equal  0.3  # Z33 
variable        sigma_c  equal  0.7  # IgG beads
variable        sigma_d  equal  0.9  # fiber center bead 
#======= average diameter over all particles =================
variable        sigma_aa  equal ${sigma_a}
variable        sigma_ab  equal 0.5*(${sigma_a}+${sigma_b})
variable        sigma_ac  equal 0.5*(${sigma_a}+${sigma_c})
variable        sigma_ad  equal 0.5*(${sigma_a}+${sigma_d})
variable        sigma_bb  equal ${sigma_b}
variable        sigma_bc  equal 0.5*(${sigma_b}+${sigma_c})
variable        sigma_bd  equal 0.5*(${sigma_b}+${sigma_d})
variable        sigma_cc  equal ${sigma_c}
variable        sigma_cd  equal 0.5*(${sigma_c}+${sigma_d})
variable        sigma_dd  equal ${sigma_d}
variable        aa_cut    equal ${cut_lj}*${sigma_aa}
variable        ab_cut    equal ${cut_lj}*${sigma_ab}
variable        ac_cut    equal ${cut_lj}*${sigma_ac}
variable        ad_cut    equal ${cut_lj}*${sigma_ad}
variable        bb_cut    equal ${cut_lj}*${sigma_bb}
variable        bc_cut    equal ${cut_lj}*${sigma_bc}
variable        bd_cut    equal ${cut_lj}*${sigma_bd}
variable        cc_cut    equal ${cut_lj}*${sigma_cc}
variable        cd_cut    equal ${cut_lj}*${sigma_cd}
variable        dd_cut    equal ${cut_lj}*${sigma_dd}
#=========== define mass  ==================================
mass        1    0.1 #filler
mass        2    0.1 #ligand
mass        3    1   #linker
mass        4    6.5 #z33
mass        5    80  #fc
mass        6    80  #fab


dielectric      ${epsilon_r}
#region void cylinder z 2 3 ${radio} -${radio} EDGE units box
#=======  specify initial velocity of atoms  =================
#velocity       all create $T 102486 mom yes rot yes dist gaussian
velocity        all create ${T_mix} 8008 dist gaussian

#========= group definition ====================================
group fiber type 1 2 3 4
group IgG   type 5 6
group Fc1   type 5
group Z33_Fc type 4 5 

group rigid  type 1 2 5 6
group nonrigid type 3 4

#===== specify parameters between atoms of type i with j====
# Partcles (1,2,3 = a; 4 = b; 5,6,7 = c; 8 = d)
# pair_coeff    i  j  ${epsilon} ${sigma}
pair_style      lj/cut ${cut_lj} 
pair_modify     shift yes mix arithmetic
pair_coeff      1  1  1.0 ${sigma_aa} ${aa_cut}
pair_coeff      2  2  1.0 ${sigma_aa} ${aa_cut}
pair_coeff      3  3  1.0 ${sigma_aa} ${aa_cut}
pair_coeff      4  4  1.0 ${sigma_bb} ${bb_cut}
pair_coeff      5  5  1.0 ${sigma_cc} ${cc_cut}
pair_coeff      6  6  1.0 ${sigma_cc} ${cc_cut}
pair_coeff      7  7  1.0 ${sigma_cc} ${cc_cut}
pair_coeff      8  8  1.0 ${sigma_dd} ${dd_cut}
pair_coeff      9  9  1.0 ${sigma_cc} ${cc_cut}
pair_coeff      10 10 1.0 ${sigma_bb} ${bb_cut}
pair_coeff      3  8  1.0 0.55  0.62

#============ specify bond / angle parameters =================================
bond_style harmonic
#              type  k   length
bond_coeff      1    300   0.1       # ligand-linker
bond_coeff      2    10    0.3       # linker-Z33 
bond_coeff      3    30    0.53      # Z33 - Fc COM distance

angle_style harmonic
angle_coeff 1 300 135 

special_bonds lj/coul 0 0 0

#  specify parameters for neighbor list 0.3 bin for units = lj, 2.0 bin for units = real or metal
neighbor       0.3 bin
comm_modify     cutoff 8

write_data      init.data  pair ij

#======= VERLET INTEGRATION WITH LANGEVIN THERMOSTAT  =======
# fix ID group-ID langevin Tstart Tstop damp seed keyword
fix             fnve rigid rigid/nve/small molecule
fix             fnve2 nonrigid nve
compute         dtemp1 fiber temp
compute_modify  dtemp1 dynamic yes
compute         dtemp2 IgG temp
compute_modify  dtemp2 dynamic yes

fix             fT1 fiber langevin ${temp}  ${temp} 10.0 123
fix             fT2 IgG   langevin ${temp}  ${temp} 100.0 16516
fix_modify      fT1 temp dtemp1 
fix_modify      fT2 temp dtemp2

#=============== mix the system first; no reaction =======================
thermo         500
thermo_style custom step pe c_dtemp1 c_dtemp2 ke evdwl epair ebond enthalpy
timestep        0.001
run             1000
write_data      EQed.data pair ij
reset_timestep  0
#=============== create bonds =============================================

molecule         mol1 step1_unreact.molecule
molecule         mol2 step1_react.molecule

fix              react1 all bond/react reset_mol_ids no &
                 react step1 all 1 0.5 0.6 mol1 mol2 step1_map prob 1 1651 


compute        bond Z33_Fc property/local batom1 batom2 btype

thermo         400
thermo_style custom step pe c_dtemp1 c_dtemp2 f_bind1[*] f_bind2[*]

dump    bonds_dump all local 400 bonds_gene.dump index c_bond[*]
dump    a all custom 20000 md0.lammpstrj   id mol type q  x  y  z
9

You can always upload files to a folder in some cloud service (Dropbox, Google Drive, MS One Drive, etc.) and provide a shareable link to that folder.

It would also be extremely helpful to create a minimum “reproducer” input that uses a very small system and has an input from which all unnecessary commands are removed. That will speed up debugging enormously. Actually, it is even a good idea for users that want to try out a new command or keyword. :wink:

10

Hi Jacob and Axel!

As Axel suggested, I uploaded all files to the folder with the link below. Also, I generate a very small system.

bond/react/forum - Google Drive.

Thank you very much for your patience and help!
Yang

11

The input deck you linked has multiple errors related to reading in the data file. Please correct these if you want anyone to look at your system.

But, here are a couple initial comments: 1) For systems that define an ‘angle style’ (but not a ‘dihedral style’), you need to include all atoms at least two bonds away from reacting atoms. 2) Fix bond/react is known to not be compatible with rigid bodies in many cases. Consider avoiding rigid bodies if you want to use fix bond/react.