(no subject)

Dear all,

I am trying to set calculations for a dna molecule by treating each nucleobase as a separate rigid body.
For this reason I use the configuration file below:

/********************************/
units real

atom_style hybrid full dipole
bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style harmonic
pair_style hybrid/overlay lj/class2 10.0 coul/cut 10.0
dielectric 1
read_data 10bp_0.topo
pair_coeff * *12 lj/class2 1.0 1.0
pair_coeff * * coul/cut

neighbor 2.0 bin

timestep 1.0

thermo_style multi
thermo 1

group rigids type 4 5 6 7 8 9 10 11 13
group no-rigids subtract all rigids

#delete_bonds rigids multi
#neigh_modify exclude molecule rigids

fix 1 no-rigids nvt temp 1.0 1.0 5.0
fix 2 rigids rigid/nvt molecule temp 1.0 1.0 5.0

run 0

/*****************************************/

So the types 4 5 6 7 8 9 10 11 13 are the types of the atoms in the nucleobase => rigids group
and the rest of the nucleoside => no-rigids group

I have to fixes. The first one responsible to integrate the no-rigids group works with no problem if it is the only fix. When I add the next fix (2) I take the problem below.

//
LAMMPS (5 May 2012)
Scanning data file …
2 = max bonds/atom
3 = max angles/atom
5 = max dihedrals/atom
1 = max impropers/atom
Reading data file …
orthogonal box = (-236.84 -241.014 -229.768) to (264.212 268.385 257.139)
1 by 1 by 1 MPI processor grid
230 atoms
148 bonds
204 angles
248 dihedrals
10 impropers
Finding 1-2 1-3 1-4 neighbors …
3 = max # of 1-2 neighbors
4 = max # of 1-3 neighbors
8 = max # of 1-4 neighbors
10 = max # of special neighbors
170 atoms in group rigids
60 atoms in group no-rigids
20 rigid bodies with 170 atoms /
which is correct since I have 20 nucleosides Also I checked and the number of atoms is correct */
Segmentation fault
/
******/

I would appreciate any suggestion about what is going wrong.

Kind regards,
Chris

Hi Chris,

could you try fix rigid with langevin to see if the seg fault occurs?
Could you upload the data file, too?

Regards,
-Trung

In all of the integrators I have the same problem. Moreover, when I use
group rigids type 13
it works. Type 13 corresponds to atoms which are not bonded to any other atom. All other atom types in rigid group correspond to atoms which are bonded within the rigid body. Also each rigid body (nucleobase) is bonded to the rest of the model (phosphate and sugar)

Many Thanks

Chris

In all of the integrators I have the same problem. Moreover, when I use
group rigids type 13
it works. Type 13 corresponds to atoms which are not bonded to any other atom. All other atom types in rigid group correspond to atoms which are bonded within the rigid body. Also each rigid body (nucleobase) is bonded to the rest of the model (phosphate and sugar)

this looks like it needs some in depth review of the affected code.
in order to resolve this properly, please produce a small test case
and post it here. people that have the experience to debug such
issues, rarely have the time to set something like this up by themselves.

thanks,
     axel.

Many thanks for the reply,

Below is a test case where a take a segmentation fault as I described before.

Configuration file:
/**********************************************************/

units real

atom_style hybrid full dipole
bond_style harmonic
angle_style harmonic
#dihedral_style opls
#improper_style harmonic
pair_style hybrid/overlay lj/class2 10.0 coul/cut 10.0
dielectric 1
read_data test.data
pair_coeff * *2 lj/class2 1.0 1.0
pair_coeff * * coul/cut

neighbor 2.0 bin

timestep 1.0

thermo_style multi
thermo 1

group rigids type 3
group no-rigids subtract all rigids

#delete_bonds rigids multi
#neigh_modify exclude molecule rigids

fix 1 no-rigids nvt temp 1.0 1.0 5.0
fix 2 rigids rigid/nvt molecule temp 1.0 1.0 5.0

run 0

/**************************************************************/

and test.data file

/**************************************************************

LAMMPS Description

         6 atoms
         4 bonds
         5 angles

         3 atom types
         1 bond types
         1 angle types

-16.840194 15.211560 xlo xhi
-11.013691 18.385058 ylo yhi
-19.768095 17.139462 zlo zhi

Masses

  1 4.0
  2 3.5
  3 0.1

Bond Coeffs

  1 1 3.48

Angle Coeffs

  1 1 120.0

Atoms

  1 1 0.00 0.00 0.00 A1 1 0.000 1.0 1.00 0.0
  2 2 1.00 1.00 1.00 A1 1 0.000 1.0 2.00 0.0
  3 2 1.50 3.00 1.00 A1 1 0.000 1.0 3.00 0.0
  4 2 1.50 0.00 1.00 A1 1 1.000 0.0 4.00 0.0
  5 3 3.00 3.00 2.00 A1 1 -1.000 0.0 5.00 0.0
  6 3 3.00 -2.00 2.00 A1 1 1.000 0.0 6.00 0.0
Bonds

  1 1 1 2
  2 1 2 3
  3 1 2 4
  4 1 3 4

Angles

  1 1 1 2 3
  2 1 1 2 4
  3 1 2 3 4
  4 1 2 4 3
  5 1 3 2 4

/****************************************************************************/

Again thanks for helping resolve this issue and for your time.

Chris

Hi Chris,

I suppose there's a bug in fix rigid involving getting atom mass.
Could you try recompiling LAMMPS with the attached fix_rigid.cpp to
see if it works?

Cheers,
-Trung

diff fix_rigid.cpp
714,715c714
< if (rmass) massone = rmass[i];
< else massone = mass[type[i]];

fix_rigid.cpp (63.5 KB)

Dear Trung,

Now it seems to work fine.

Thanks,
Chris

Good catch Trung - I'll post this as a patch.
Must have missed these 2 spots for checking
for per-type vs per-atom mass when extending
fix rigid to work with both types of particles.

Steve