conjugate gradient minimization

Moumita_Maiti · August 13, 2015, 11:40am

Hi,

I am minimizing soft repulsive harmonic potential using conjugate gradient method in lammps. The line minimization I use is quadratic. I have added a small part to calculate force and energy for this potential as lammps have no source code with this potential. I am trying with four particles in a box of length 2, so the minimized energy would be zero.

I face following problem:

My code gives minimized energy zero but lammps stop before minimizing. It is very sensitive to dmax value mentioned in input file for line minimization.

I figured out one reason for that, I see I can not mention dmax value for me greater than 0.3, whereas zero energy is found if dmax is minimum 0.8. If I use more than 0.3, it calculates wrong energy which I do not understand.

Please help me!

Regards,

Moumita

akohlmey · August 13, 2015, 2:27pm

Hi,

I am minimizing soft repulsive harmonic potential using conjugate gradient
method in lammps. The line minimization I use is quadratic. I have added a
small part to calculate force and energy for this potential as lammps have
no source code with this potential. I am trying with four particles in a
box of length 2, so the minimized energy would be zero.

I face following problem:

My code gives minimized energy zero but lammps stop before minimizing. It is
very sensitive to dmax value mentioned in input file for line minimization.

I figured out one reason for that, I see I can not mention dmax value for
me greater than 0.3, whereas zero energy is found if dmax is minimum 0.8. If
I use more than 0.3, it calculates wrong energy which I do not understand.

Please help me!

it is near impossible to give advice based on such vague and indirect
descriptions and specifically no example input deck.

axel.

Moumita_Maiti · August 13, 2015, 2:33pm

Well, my input script is as follows:

akohlmey · August 13, 2015, 2:35pm

Well, my input script is as follows:

#####################
dimension 3
newton on
comm_modify vel yes
read_data initcoord.in
group g1 type 1
mass * 1.0
pair_style dpd 0.000000 1.000000 478888
pair_coeff 1 1 1.000000 1.000000 1.000000
neighbor 1.750000 bin
neigh_modify every 1 delay 0 check no
thermo 10
velocity all set 0. 0. 0. units box
min_style cg
min_modify dmax 0.2 line quadratic
minimize 0.0000000001 10000 10000
dump mydumpCoord all custom 1 dump.10000.10 id type x y z
dump_modify mydumpCoord format "%d d .16g \.16g .16g"
thermo_style custom step press temp pe
fix nve1 all nve
run 1

Input configuration is:

#check...
4 atoms
1 atom types
0.0 2.0 xlo xhi
0.0 2.0 ylo yhi
0.0 2.0 zlo zhi
0.0 0.0 0.0 xy xz yz
Atoms

1 1 0.0 0.0 0.0
2 1 0.5 0.2 0.0
3 1 0.5 0.0 0.2
4 1 0.0 0.5 0.2

I have modified pair_dpd.cpp source code for repulsive harmonic potential.

nobody can debug problems on modified code.
how do you know that the problem isn't due to a mistake in your
modifications?
i guess your are on your own on this one.

axel.

Moumita_Maiti · August 13, 2015, 2:44pm

Thanks! I am trying that. I thought to share the piece of code I modified which is pasted below. I checked each term by term energy and force and it is doing correct.

What I found is that dmax sets maximum distance a particle can move in the line minimization. If a particle needs to move larger than the given dmax value in order to find minimum in line-minimization, it can not, rather it returns with a higher value. I could check this as I have my own minimization code and comparing results of both from lammps and my code. I need lammps to work because I need to run for very large system size. Anyway thanks.

Regards,
Moumita
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
vxtmp = v[i][0];
vytmp = v[i][1];
vztmp = v[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];

for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_dpd = special_lj[sbmask(j)];
j &= NEIGHMASK;

delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delz*delz;
jtype = type[j];

if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
rinv = 1.0/r;
delvx = vxtmp - v[j][0];
delvy = vytmp - v[j][1];
delvz = vztmp - v[j][2];
dot = delxdelvx + delydelvy + delz*delvz;
wd = 1.0 - r/cut[itype][jtype];
randnum = random->gaussian();

f[i][0] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delzfpair;
printf(“i=%d %f %f %f\n”, i, f[i][0],f[i][1],f[i][2]);
if (newton_pair || j < nlocal) {
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
f[j][2] -= delzfpair;
printf(“j=%d %f %f %f\n”, j, f[j][0],f[j][1],f[j][2]);
}

if (eflag) {
// unshifted eng of conservative term:
// evdwl = -a0[itype][jtype]r * (1.0-0.5r/cut[itype][jtype]);
// eng shifted to 0.0 at cutoff
evdwl = wd*wd;
//evdwl *= factor_dpd;
}

if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
}

Moumita_Maiti · August 13, 2015, 2:56pm

In case somebody follows this message, here is the correct code. In last message I deleted two lines where fpair is calculated.

Regards,
Moumita

for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
vxtmp = v[i][0];
vytmp = v[i][1];
vztmp = v[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];

for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_dpd = special_lj[sbmask(j)];
j &= NEIGHMASK;

delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delz*delz;
jtype = type[j];

if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
rinv = 1.0/r;
delvx = vxtmp - v[j][0];
delvy = vytmp - v[j][1];
delvz = vztmp - v[j][2];
dot = delxdelvx + delydelvy + delz*delvz;
wd = 1.0 - r/cut[itype][jtype];
randnum = random->gaussian();

fpair = 2.0*wd;
fpair = 1./cut[itype][jtype];fpair=fpairrinv;

f[i][0] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delzfpair;
printf(“i=%d %f %f %f\n”, i, f[i][0],f[i][1],f[i][2]);
if (newton_pair || j < nlocal) {
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
f[j][2] -= delzfpair;
printf(“j=%d %f %f %f\n”, j, f[j][0],f[j][1],f[j][2]);
}

if (eflag) {
// unshifted eng of conservative term:
// evdwl = -a0[itype][jtype]r * (1.0-0.5r/cut[itype][jtype]);
// eng shifted to 0.0 at cutoff
evdwl = wd*wd;
//evdwl *= factor_dpd;
}

if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
}

if (vflag_fdotr) virial_fdotr_compute();

akohlmey · August 13, 2015, 3:15pm

Thanks! I am trying that. I thought to share the piece of code I modified
which is pasted below. I checked each term by term energy and force and it
is doing correct.

it is still a bad example since, a) it is always better to write a new
style rather than hacking an existing one (the new one can still be based
on the old), otherwise input are *far* to confusing and there is a very
high risk of accidentally using the wrong executable without the
modification or to accidentally overwrite the modification. and b) if you
care so much about performance, you should have used a template style that
does not require communication of velocities (let alone then do computation
with that data without using the result of that. since your computation is
extremely "cheap" such "useless" computations can significantly impact
performance).

What I found is that dmax sets maximum distance a particle can move in the
line minimization. If a particle needs to move larger than the given dmax
value in order to find minimum in line-minimization, it can not, rather it
returns with a higher value. I could check this as I have my own
minimization code and comparing results of both from lammps and my code. I
need lammps to work because I need to run for very large system size.
Anyway thanks.

what you say doesn't make sense. why should you need to set such a large
dmax? the whole point of limiting the maximum displacement is to avoid bad
situations when particles overlap too much. at best, it will take some
additional iterations to get to the minimum.

in any case, a quick read through the documentation for min_modify reveals
that for your system and requirements, you should use the forcezero
linesearch instead of quadratic. if i use this with a purely repulsive
lj/cut style using a 1 sigma cutoff, it finds the minimum quickly and
reliably.

also, your input doesn't throws an error on the minimize command. do you
have more modifications in you version of LAMMPS that you didn't tell us
about?

axel.

Moumita_Maiti · August 13, 2015, 3:35pm

Thanks! I am trying that. I thought to share the piece of code I modified
which is pasted below. I checked each term by term energy and force and it
is doing correct.

it is still a bad example since, a) it is always better to write a new
style rather than hacking an existing one (the new one can still be based
on the old), otherwise input are *far* to confusing and there is a very
high risk of accidentally using the wrong executable without the
modification or to accidentally overwrite the modification. and b) if you
care so much about performance, you should have used a template style that
does not require communication of velocities (let alone then do computation
with that data without using the result of that. since your computation is
extremely "cheap" such "useless" computations can significantly impact
performance).

Thanks, yes I hacked as it was easy for me. I like point b), I will
surely try this.

What I found is that dmax sets maximum distance a particle can move in
the line minimization. If a particle needs to move larger than the given
dmax value in order to find minimum in line-minimization, it can not,
rather it returns with a higher value. I could check this as I have my own
minimization code and comparing results of both from lammps and my code. I
need lammps to work because I need to run for very large system size.
Anyway thanks.

what you say doesn't make sense. why should you need to set such a large
dmax? the whole point of limiting the maximum displacement is to avoid bad
situations when particles overlap too much. at best, it will take some
additional iterations to get to the minimum.

in any case, a quick read through the documentation for min_modify reveals
that for your system and requirements, you should use the forcezero
linesearch instead of quadratic. if i use this with a purely repulsive
lj/cut style using a 1 sigma cutoff, it finds the minimum quickly and
reliably.

also, your input doesn't throws an error on the minimize command. do you
have more modifications in you version of LAMMPS that you didn't tell us
about?

No additional modification. Seems if energy is zero in line-minimization,

still it does not stop. So, I did not get any error but higher minimized
energy. When I added this modification in the line_minimization code that
if energy is zero, stop and go back to conjugate gradient code, it gave me
correct output configuration with zero energy for larger dmax 0.3.

Now, coming to the dmax value, you are correct it takes longer time with
low value of dmax, and minimizes. But it does not end in same
configuration, the problem is there. I calculate number of contacts at
lower density and this is very sensitive to dmax. Where contact number
should be zero, I get around 3, as it kept overlap of the order of 10^-12
or so. Lets see. But thanks a lot for so much input.

Regards,
Moumita

Moumita_Maiti · August 14, 2015, 7:59am

As you suggested Axel, forcezero linesearch instead of quadratic works quickly and reliably for purely repulsive potential. It does, Thanks!

Regards,
Moumita

sjplimp · August 14, 2015, 3:07pm

If you are doing minimization with pair DPD, you probably

need to turn off the random force term. I don’t think the

minimizer will work if that is added in every time forces

are evaluated.

Steve