[lammps-users] energy minimization + constraint


I have a question concerning the energy minimization (using cg) using constraints.
I figured out that the results of energy minimization are very meaningful, when there are no constraints at all.

My simulation box is 3d with (s) as boundary conditions in x,y-direction and § periodic boundary condition in z-direction
I apply forces in y-dir. at the top and the bottom atoms of the system by using the fix setforce.

If I apply a zero force, the results are very meaningful, basically the concerned atoms are fixed.
If I apply a force, which is unequal to zero, then non-reasonable things happen (I calculate, e.g., the virial stress at the centre of the simulation box after I brought the system to equilibrium by energy minimization).
If I apply forces from f_1=10 to f_10=100, the overall trend of the response of stress is correct (the stresses increase), but it is by no means a linear increase. It might happen that f_3 gives nearly the same results as f_4 and then there is a big jump to f_5.

My questions:

  1. Is there anybody out there, who ever did an energy minimization with force constraints, which are unequal to zero?
  2. Can anybody point me out a paper or description of the conjugate gradient method as energy minimization (I know about CG and solving algebra systems, but the algorithm in LAMMPS should be a bit different)

Thank you.

Using non-zero forces with fix setforce should work. I've played
with it for a few simple cases. I would try visualizing what
your system is doing as you ramp up the force. It may be
that you are introducing other defects or changes to your
system that are non-linear. How are you measuring the
virial at the center of the box?

The src/min.cpp file lists this citation which is very good.
"Conjugate Gradient Method Without the Agonizing Pain" by
            JR Shewchuk, http://www-2.cs.cmu.edu/~jrs/jrspapers.html#cg


I apply to all upper and lower atoms a force and calculate by using fix ave/spatial a qualitative measure for the stress (with stress/atom) at the centre of the system (dimension is [bar*Angstrom^3]). This can be transformed to a virial stress, however, just by itself it gives an idea for the stress for comparison.
There are no cracks or flaws in the system.

If I apply different loads in y-direction, I get following behaviour (please, watch sigma_yy)

no_atoms sig_xx sig_yy sig_zz (pos top atoms)
0 force
108 -10953.6 -23048 -14525.5 72.3
0.000025 force
108 -12884.7 -22349.5 -15016 72.3048 +0.0066%
0.000050 force
108 -13328.7 -22220.2 -15155.9 72.3094 +0.0130%
0.000075 force
108 -12735.1 -20460.6 -14088.8 72.3134 +0.0185%
0.000100 force
108 -15705.1 -21456.1 -15829.9 72.3168 +0.0232%

(no_atoms is the number of atoms at the centre out of which the stress is calculated -> stays the same)
(pos_top_atoms gives the position in y-direction of the top atoms -> behaves meaningfully)
all shear stresses are zero, as they should be

This behaviour of the normal stress in force direction is strange to me and I cannot explain it (the strain behaves nearly linear, the stress not at all).

I attach my lammps program for reference how I calculate the stress, it is too hard to explain in words,
Thank you,

I was wrong. Fix setforce is not setup to be used
with minimization when the applied force is non-zero.
The doc page should say this explicitly. The reason
is that there is no energy term associated with
that force, which would need to be done by integrating
over the distance the force acts.

I think you want to use fix addforce which does
this, as the energy explanation in the doc page