[lammps-users] problem with "compute dipole/chunk"

Dear lammps users,

I want to compute the dipole moment of a spherical region in the bulk liquid water. I test two routes to obtain the dipole moment:

1 . Define a spherical region, and make a dynamic group to put all atoms in this spherical region volume into the group, then use “variable Mx atom x*q” and “compute reduce sum” to get dipole moment.

  1. Use “chunk/atom” to define a spherical region, and compute the dipole moment by “compute dipole/chunk”.

Please find more details in the test.in file.

I print the results every 100 steps by:

fix test1 all ave/time 100 1 100 c_dipole1[*] file tmp1.out mode vector

fix test2 all ave/time 100 1 100 c_dipole[*] file tmp2.out mode vector

My problem is those two routes can not give the same values of the dipole moment. In addition, as it is mentioned in the manual about “compute dipole/chunk” command, one may reset the image flags (e.g. to 0) before invoking this compute by using the “set image” command. However, there is no problem if I want to compute the total dipole moment of the whole bulk system, but it gives weird results for a spherical region.

I am looking forward to your reply

Best wishes,

Wei

tmp2.out (152 Bytes)

tmp1.out (200 Bytes)

test.in (1.6 KB)

system.in.init (312 Bytes)

system.data (655 KB)

system.in.settings (222 KB)

you are likely not taking into account what it is exactly that compute dipole/chunk computes. specifically, there is:

For chunks with a net charge the resulting dipole is made position independent by subtracting the position vector of the center of mass or geometric center times the net charge from the computed dipole vector.

so your computation of: Sum_n q_ir_i is not what compute dipole/chunk is computing which is instead (by default): Sum_n q_ir_i - Sum_n q_i * (Sum_n m_ir_i)/Sum_n m_i
or Sum_n q_i
r_i - Sum_n q_i * Sum_nr_i (when using the “geom” option).

when cutting a spherical region of water molecules, you will likely have a net charge. with water molecules, the cleaner solution in either case would be to make the group selection or chunk assignment such that whole molecules are selected where then the ambiguity and position dependence of Sum_n q_i*r_i due to a net charge is avoided by construction.

Axel.

Dear Axel,

Thank you very much for your explanation.

I still have one more question about the command “set image”. The image flag is reset before the “compute dipole/chunk” command, the value of the dipole moment is much greater than that without "set image"command. While there is no problem if I want to compute the total dipole moment of the whole system.

I am looking forward to your reply.

Best wishes,

Wei

you must not blindly just reset all image flags to zero, especially not for molecules. in cases where a molecule stretches across a periodic boundary, you would break it into pieces as normally part of the molecule would have different image flags than other parts.

for neutral molecules in general, it would be best (for total dipole moment and per chunk/group/region dipole moment) to do the selection on just one reference atom of the molecule based on wrapped coordinates (i.e. without consideration of image flags) and then selecting the whole molecule as then the computation of the dipole moment should use unwrapped coordinates and will be position independent and more meaningful.

if there needs to be a reset of image flags, one might try writing out a data file and then processing the image flags in there with a custom script. or do it on a dump file and the first read the data file and then update the positions (and image flags) with read_data.

axel.

Axel

sorry, but what you are telling us is getting more inconsistent with every e-mail and thus leading nowhere.
you need to start over and provide a suitable, small but complete input deck and a specific explanation of where you see a problem and describe in detail how to reproduce it.
i don’t see any problems with my examples and cannot tell whether you are doing things correctly or if the situation is exactly as you claim it is.

axel.

Dear Axel,

Actually, the volume of the spherical region is smaller than the whole simulation box. The molecules in this sphere must be complete, none of them across the periodic boundary. So they all have the same image flags, “set image” should not affect the result.

And, the dipole moment should increase with increasing volume(M=P*V, M: dipole moment, P:polarization, V:volume). However, as I change the radius of the spherical region, the “compute dipole/moment” always give the same result, even the number after the decimal point is the same.

R=5.9 A,

Row c_dipole1[1] c_dipole1[2] c_dipole1[3] c_dipole1[4]

1000000 1
1 41.3806 -13.1636 46.0215 63.2741

R=11.8A

Row c_dipole1[1] c_dipole1[2] c_dipole1[3] c_dipole1[4]

1000000 1
1 41.3806 -13.1636 46.0215 63.2741

R=17.7A

Row c_dipole1[1] c_dipole1[2] c_dipole1[3] c_dipole1[4]

1000000 1
1 41.3806 -13.1636 46.0215 63.2741

Best wishes,

Wei

Dear Axel,

The previous email is too big, I compress it and send it again.

I want to compute the dipole moment of a spherical region in the bulk liquid water. I use “chunk/atom” to define a spherical region,

variable Cx equal 28.0 # the center of the sample in the x direction

variable Cy equal 28.0 # the center of the sample in the y direction

variable Cz equal 28.0 # the center of the sample in the z direction

variable R equal 19.64157993 # the radius of the spherical vessel

variable rc equal 0.5*v_R

compute mychunk all chunk/atom bin/sphere {Cx} {Cy} {Cz} 0.0 {rc} 1 nchunk every ids every

and compute the dipole moment by “compute dipole/chunk”.

compute dipole1 all dipole/chunk mychunk

I print the results every 100 steps by:

fix test1 all ave/time 100 1 100 c_dipole1[*] file tmp1.out mode vector

Please find more details in the test.in file.

I find two major problems:

  1. The image flag is reset before the “compute dipole/chunk” command, the value of the dipole moment is much greater than that without "set image"command. While there is no problem if I want to compute the total dipole moment of the whole system. The volume of the spherical region is smaller than the whole simulation box. The molecules in this sphere must be complete, none of them across the periodic boundary. So they all have the same image flags, “set image” should not affect the result.

  2. The dipole moment should increase with increasing volume(M=P*V, M: dipole moment, P:polarization, V:volume). However, as I change the radius of the spherical region ${rc}, the “compute dipole/moment” always give the same result, even the number after the decimal point is the same.

rc =5.9 A,

Row c_dipole1[1] c_dipole1[2] c_dipole1[3] c_dipole1[4]

1000000 1
1 41.3806 -13.1636 46.0215 63.2741

rc**=11.8A**

Row c_dipole1[1] c_dipole1[2] c_dipole1[3] c_dipole1[4]

1000000 1
1 41.3806 -13.1636 46.0215 63.2741

rc**=17.7** A

Row c_dipole1[1] c_dipole1[2] c_dipole1[3] c_dipole1[4]

1000000 1
1 41.3806 -13.1636 46.0215 63.2741

Best wishes,

Wei

inscript.tar.xz (90.7 KB)