I am trying to make rigid a CO2 molecule with fix shake. but when i use the : fix 1 all shake 1e-4 500 0 b 1 a 1 i get the error ERROR: Out of range atoms - cannot compute PPPM.

I cannot find any mistake in my command.I am trying to change it with fix 1 all shake 1e-4 500 0 m mVALUE a 1 but in the mValue do i have to put the C mass? Do i make any other mistake Because i have simulate the same prob without kspace and a pair_style lj/cut/coul/cut and there is no problem.

Regards

Dimitris

dimension 3

boundary p p p

units real

neighbor 2 bin

neigh_modify delay 2 every 1

atom_style full

bond_style harmonic

angle_style harmonic

#pair_style lj/cut/coul/cut 8 10

pair_style lj/cut/coul/long 8 10

pair_modify mix arithmetic

kspace_style pppm 1e-4

read_restart CO2data.minim-Harris.restart

pair_coeff 1 1 0.159983896019 3.033 #OO epsilon(kcal/mol) sigma(A)

pair_coeff 2 2 0.0558980835347 2.757 #CC epsilon-sigma

bond_coeff 1 10000 1.149

angle_coeff 1 10000 180.0

fix 1 all shake 1e-4 500 10 b 1 a 1

fix 2 all npt temp 300 300 100.0 iso 0.98 0.98 1000.0

#equilibrate

thermo 300

thermo_style multi

timestep 0.5

run 200000

write_restart co2.equil.restart

#production run

thermo 500

log produc.log

thermo_style custom etotal ke temp pe ebond eangle ecoul press vol

thermo_modify line multi

timestep 1.0

run 1000000

Does fix SHAKE tell you that it applied the correct constraints to one

or more CO2 molecules?

Steve

Yes it does. Moreover I change it to fix 1 all shake 1e-04 500 0 t 1 2

and it works fine and i still have the molecule rigid based on the manual.But why this is happening?

I don't understand what you are asking. You seem to be saying

that you specify the SHAKE command in two different ways,

and one of them works fine, but the other one blows up. Yet

in both cases, the SHAKE output is telling you that it has chosen

to constrain the same number of bonds/angles, etc, and if

you print the SHAKE stats out, they are the same?

If that is the case, it seems unlikely the 2 simulatoins will be

different. The PPPM error you are getting typically means

the simulation has blown up, and atoms are flying far apart.

Steve

steve,

does the fix shake in LAMMPS have a special case treatment

for linear molecules? i don't see it and thus i think it cannot work.

the way i understand the code, it converts the constraint atom selection

into two bond and one angle constraint. the angle constraint at 180

degrees cannot be resolved, as the forces are degenerate. which in

turn would lead to Inf or NaN in the shake forces and thus kill pppm

as a collateral victim.

cheers,

axel.

Possibly this is an issue. But I still don't understand what is different

about the 2 cases. The args for the fix shake command simply

enable the code to choose what to constrain. If the 2 commands

produce the same list of constraints, then the subsequent behavior

should be identical. Dimitrios - please post the output of

the 2 runs, including what SHAKE says it is constraining, and the

first few steps of thermo output and turn on SHAKE stats. I want

to see if anything is different.

Steve

WITH fix 1 all shake 1e-4 500 10 t 1 2

LAMMPS (14 Apr 2010)

Reading restart file ...

orthogonal box = (-0.00610646 -0.00610646 -0.00610646) to (347.782 347.782 347.782)

2 by 2 by 2 processor grid

3063 atoms

2042 bonds

1021 angles

Finding 1-2 1-3 1-4 neighbors ...

2 = max # of 1-2 neighbors

1 = max # of 1-3 neighbors

1 = max # of 1-4 neighbors

2 = max # of special neighbors

Finding SHAKE clusters ...

0 = # of size 2 clusters

1021 = # of size 3 clusters

0 = # of size 4 clusters

0 = # of frozen angles

PPPM initialization ...

G vector = 0.110536

grid = 40 40 40

stencil order = 5

RMS precision = 8.82157e-05

brick FFT buffer size/proc = 15625 8000 5625

Setting up run ...

SHAKE stats (type/ave/delta) on step 95

1 1.14908 0.0298712

Memory usage per processor = 4.92145 Mbytes

---------------- Step 95 ----- CPU = 0.0000 (sec) ----------------

TotEng = 3365.0791 KinEng = 1681.8856 Temp = 236.9419

PotEng = 1683.1935 E_bond = 0.0000 E_angle = 1699.0799

E_dihed = 0.0000 E_impro = 0.0000 E_vdwl = -14.7502

E_coul = 13446.8464 E_long = -13447.9826 Press = -28.7446

Volume = 42067117.3115

SHAKE stats (type/ave/delta) on step 100

1 1.149 7.39937e-10

SHAKE stats (type/ave/delta) on step 110

1 1.149 2.09058e-09

SHAKE stats (type/ave/delta) on step 120

1 1.149 5.67342e-09

SHAKE stats (type/ave/delta) on step 130

1 1.149 2.1964e-09

SHAKE stats (type/ave/delta) on step 140

1 1.149 4.76849e-09

SHAKE stats (type/ave/delta) on step 150

1 1.149 1.13929e-08

SHAKE stats (type/ave/delta) on step 160

1 1.149 3.6592e-09

SHAKE stats (type/ave/delta) on step 170

1 1.149 8.37158e-09

SHAKE stats (type/ave/delta) on step 180

1 1.149 1.8152e-08

SHAKE stats (type/ave/delta) on step 190

1 1.149 5.20014e-09

SHAKE stats (type/ave/delta) on step 200

1 1.149 1.22928e-08

SHAKE stats (type/ave/delta) on step 210

1 1.149 2.29537e-08

SHAKE stats (type/ave/delta) on step 220

1 1.149 5.78698e-09

SHAKE stats (type/ave/delta) on step 230

1 1.149 1.40289e-08

SHAKE stats (type/ave/delta) on step 240

1 1.149 2.22527e-08

SHAKE stats (type/ave/delta) on step 250

1 1.149 4.84604e-09

SHAKE stats (type/ave/delta) on step 260

1 1.149 1.28006e-08

SHAKE stats (type/ave/delta) on step 270

1 1.149 1.76721e-08

SHAKE stats (type/ave/delta) on step 280

1 1.149 3.47799e-09

SHAKE stats (type/ave/delta) on step 290

1 1.149 1.04366e-08

SHAKE stats (type/ave/delta) on step 300

1 1.149 1.29648e-08

---------------- Step 300 ----- CPU = 3.7310 (sec) ----------------

TotEng = 2935.0574 KinEng = 2573.3591 Temp = 362.5316

PotEng = 361.6983 E_bond = 0.0000 E_angle = 391.0989

E_dihed = 0.0000 E_impro = 0.0000 E_vdwl = -21.0774

E_coul = 13439.6631 E_long = -13447.9863 Press = -0.4479

Volume = 42034894.1645

fix 1 all shake 1e-4 500 10 b 1 a 1

Reading restart file ...

orthogonal box = (-0.00610646 -0.00610646 -0.00610646) to (347.782 347.782 347.782)

2 by 2 by 2 processor grid

3063 atoms

2042 bonds

1021 angles

Finding 1-2 1-3 1-4 neighbors ...

2 = max # of 1-2 neighbors

1 = max # of 1-3 neighbors

1 = max # of 1-4 neighbors

2 = max # of special neighbors

Finding SHAKE clusters ...

0 = # of size 2 clusters

0 = # of size 3 clusters

0 = # of size 4 clusters

1021 = # of frozen angles

PPPM initialization ...

G vector = 0.110536

grid = 40 40 40

stencil order = 5

RMS precision = 8.82157e-05

brick FFT buffer size/proc = 15625 8000 5625

Setting up run ...

SHAKE stats (type/ave/delta) on step 95

1 1.14908 0.0298712

1 179.333 1.61838

Memory usage per processor = 4.76886 Mbytes

---------------- Step 95 ----- CPU = 0.0000 (sec) ----------------

TotEng = 1665.9992 KinEng = 1681.8856 Temp = 276.4516

PotEng = -15.8864 E_bond = 0.0000 E_angle = 0.0000

E_dihed = 0.0000 E_impro = 0.0000 E_vdwl = -14.7502

E_coul = 13446.8464 E_long = -13447.9826 Press = nan

Volume = 42067117.3115

ERROR: Out of range atoms - cannot compute PPPM

using the flags 't' and 'm' with fix shake only constrains bonds,

not angles. the documentation is very clear about that.

what i am saying is, that if you leave away the 'a' then you have the same

as you get with 'm' and 't'. otherwise you are comparing apples and oranges.

as it wrote to steve, you cannot constrain

a 180 degree angle directly.

it is matter of elementary geometry.

Finding SHAKE clusters ...

0 = # of size 2 clusters

1021 = # of size 3 clusters

0 = # of size 4 clusters

0 = # of frozen angles

Finding SHAKE clusters ...

0 = # of size 2 clusters

0 = # of size 3 clusters

0 = # of size 4 clusters

1021 = # of frozen angles

These 2 ouptuts from the SHAKE setup are not

the same, which is why I asked my original question.

The first is constraining the 1021 water mols by

bond only (flexible angle), the latter is constraining the bonds

and angles. So they are different models. It may

be, as Axel said, that the code will have a problem with

constraining 180 degree angles. I need to check that.

But regardelss, they are different models, and you should

not expect identical behavior.

Steve