Dear users and developers,
I run some simulations in the npt ensemble using the Gay-Berne potential and setting the pressure to 2.3. As I see in the output (I guess the units should be the same as in the input file) the equilibrium pressure is not 2.3 but about 2.0. I cannot find the reason for such an outcome. I tried first with an old version (lammps-20Apr12) then updated to the last version (lammps-5Mar13) with exactly the same results. One input file and the corresponding shortened output file is below. I tried different n. of molecules (from 1728 to 4096) and again the same result. If I set a lower pression in the fix command, e.g. 2.0, I get an output pressure scaled by the same factor, that is about 1.74 as equilibrium pressure. The example below is quite a crude equilibration starting from a hand made box, but I also have much longer equilibrations at several temperatures and the pressure is always different from wha I set in the input (that is lower by about a factor 1.15). In contrast, the temperature equilibrates to the correct input value set in the fix npt command.
INPUT:
units lj
atom_style ellipsoid
dimension 3
lattice sc 0.015
region box block 0 12 0 12 0 12
create_box 1 box
create_atoms 1 box
set type 1 mass 1.0
set type 1 shape 4.4 1 1
set group all quat 1 0 0 0
velocity all create 2.4 87287 loop geom
pair_style gayberne 1.0 1.0 1.0 5.5
pair_coeff 1 1 1.0 1.0 1 1 0.05 0 0 0
neighbor 0.8 bin
thermo_style one
thermo 1000
timestep 0.0015
compute q all property/atom quatw quati quatj quatk
dump 1 all custom 1000 trj.GB &
id type x y z c_q[1] c_q[2] c_q[3] c_q[4]
fix 1 all npt/asphere temp 1.45 1.45 0.1 iso 2.3 2.3 10
run 100000
OUTPUT:
LAMMPS (5 Mar 2013)
reset 1 OpenMP thread(s) per MPI task
using OpenMP capable neighbor list subroutines
using double precision OpenMP force kernels
Lattice spacing in x,y,z = 4.0548 4.0548 4.0548
Created orthogonal box = (0 0 0) to (48.6576 48.6576 48.6576)
Created 1728 atoms
Setting atom values ...
Last active /omp style is pair_style gayberne/omp
Setting up run ...
Memory usage per processor = 4.14479 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 2.4 763.35192 0 766.94983 295.79037 115200
10000 1.5790717 -0.47342846 0 1.8938084 0.43089139 19253.716
20000 1.418367 -0.79746376 0 1.3288554 1.904492 10739.734
30000 1.4532934 -0.97907126 0 1.1996073 2.1531219 10056.522
40000 1.480031 -1.1110891 0 1.1076726 2.0013327 9925.4033
50000 1.4095266 -1.1819959 0 0.93107044 2.0898041 9803.0054
60000 1.4444888 -1.1974778 0 0.96800155 2.1535823 9787.8321
70000 1.4625576 -1.2944814 0 0.89808545 2.0691284 9733.2158
80000 1.4285046 -1.521947 0 0.61956992 1.8870487 9511.9039
90000 1.4691011 -1.6454068 0 0.55696963 1.9717258 9384.4816
100000 1.4582348 -1.643293 0 0.54279341 2.1162688 9296.8185
Loop time of 10821.3 on 1 procs for 100000 steps with 1728 atoms
Pair time (\) = 10739\.6 \(99\.2448\)
Neigh time \() = 11.312 (0.104535)
Comm time (\) = 17\.3263 \(0\.160113\)
Outpt time \() = 0.511476 (0.00472657)
Other time (%) = 52.5694 (0.485795)
Nlocal: 1728 ave 1728 max 1728 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 5313 ave 5313 max 5313 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 167403 ave 167403 max 167403 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 167403
Ave neighs/atom = 96.8767
Neighbor list builds = 1881
Dangerous builds = 3
System init for write_restart ...
Last active /omp style is pair_style gayberne/omp
Mike may want to comment on this as he has done
some modeling of GBerne systems.
Steve
A lot has changed since I did production runs in LAMMPS with asphere and NPT, but I would guess that your problem almost certainly has to do with the temperature computes used in LAMMPS. The pressure is a function of temperature and different computes can be used in LAMMPS for thermodynamic output and fixes. For example, npt/asphere creates a temp/asphere compute and uses this for pressure computation. Your thermodynamic output is probably using a pressure compute that does not include angular terms in the temperature. To verify, for example:
compute rot all temp/asphere
group spheroid type 1
variable dof equal count(spheroid)*3+3 # For BIAXIAL ellipsoids, for uniaxial, 2 degrees of freedom per ellipsoid
compute_modify rot extra ${dof}
compute my_pressure all pressure rot
Output of c_my_pressure in thermo_style custom should match your target pressure. This assumes compute_modify was used to add extra degrees of freedom to the fix for npt/asphere. I think that you might also be able to use c_1_press to output the pressure used by a fix directly (in this case with fix id "1").
NPT with asphere in LAMMPS is tricky. Be careful about the temperature computes, the associated temperature computes in pressure computes, adding the correct rotational degrees of freedom to the temperature computes, scalar pressures and nematic/smetic phases with orientational order, etc.
- Mike
Dear Mike,
thanks a lot! Indeed that was the mistake. It works fine now, both temperature and pressure fluctuate around the value set in the input.
Thank you very much everybody, I really appreciate your help,
Giacomo
I'm a bit confused by this:
compute rot all temp/asphere
group spheroid type 1
variable dof equal count(spheroid)*3+3 # For BIAXIAL ellipsoids,
for uniaxial, 2 degrees of freedom per ellipsoid
compute_modify rot extra ${dof}
compute my_pressure all pressure rot
as to why the compute_modify is needed. From the
compute temp/asphere doc page:
Only finite-size particles (aspherical or spherical) can be included
in the group. For 3d finite-size particles, each has 6 degrees of
freedom (3 translational, 3 rotational). For 2d finite-size particles,
each has 3 degrees of freedom (2 translational, 1 rotational).
IMPORTANT NOTE: This choice for degrees of freedom (dof) assumes that
all finite-size aspherical or spherical particles in your model will
freely rotate, sampling all their rotational dof. It is possible to
use a combination of interaction potentials and fixes that induce no
torque or otherwise constrain some of all of your particles so that
this is not the case. Then there are less dof and you should use the
compute_modify extra command to adjust the dof accordingly.
If you use fix npt/asphere, it will create a compute temp/asphere
and assign 6 dof per ellipsoid. Everything should just work, subject
to the NOTE caveat above. Is the total dof, and thus the
temperature, not correct for your system? I.e. why do you
need to use compute_modify extra?
Mike is correct that if you just use thermo_style temp or press to
print the temperature
and pressure, they will not be the values that fix npt is using for aspherical
particles. By default the thermo output uses compute temp and compute pressure
which are for point particles.
Steve
Sorry, when I originally wrote it, dof was handled automatically (since only ellipsoids were handled) and I didn't have to deal with this stuff. Is this correct?:
For biaxial ellipsoids
compute_modify rot extra 0
For uniaxial ellipsoids in 3d
variable dof equal count(spheroid)+1
compute_modify rot extra ${dof}
For uniaxial ellipsoids in 2d
compute_modify rot extra 0
If so, I think that examples/in.ellipse.gayberne needs to be updated. Thanks. - Mike
Dear Steve and Mike,
thanks again for your help. I tried both solutions suggested by Mike and I get identical results, as far as I print "my_pressure" also the pressure is correct regardless the way the extra dof are defined. Both temperature and pressure seems to converge to the correct values. I'll see tomorrow after a longer run whether this is really true.
Best wishes,
Giacomo
log.lammps_1 (16 KB)
log.lammps_2 (12 KB)
I don't think of it in terms of uniaxial vs biaxial, but rather whether
the particles are fully rotating or not. For 3d particles, if their omega
is a general 3-vector, then there are 6 dof per particle. If the force
field does not induce a full omega (i.e. 1,2, or 3 components are missing),
then you need to subtract additional dof from 6. So are you saying
that uniaxial GB particles only have 2 components of omega, and thus
5 dof per particle?
Steve
I apologize for confusing everyone including myself here.
Giacomo, for uniaxial ellipsoids, you need to subtract a degree of freedom for every uniaxial particle for every temp/asphere compute:
compute_modify rot extra ${dof}
where rot is the temperature compute, and dof equals the number of uniaxial ellipsoids. You should see different results when using this versus when not. There is a temp/asphere compute that is used by the npt/asphere class. It has an id of FIX_temp where FIX is the id of the npt/asphere fix. So to alter the degrees of freedom used in the temperature computation for the thermostat, for example:
fix 1 all npt/asphere temp 2.4 2.4 0.5 iso 0.0 8.0 0.5
compute_modify 1_temp extra ${dof}
If you are doing this and still see identical results, I think that something is wrong.
Steve,
The answer to your question is "Yes". Just to be specific, typically (I have never seen otherwise), uniaxial ellipsoids in GB are parametrized so that there will never be torque around the optical axis. This is required in the original Gay-Berne model that only handled uniaxial ellipsoids. For the berardi biaxial generalization in LAMMPS this will be the case when exactly two semiaxes have the same length and the corresponding relative well depths are equal.
- Mike
Dear Mike and Steve,
yes, with this additional compute_modify 1_temp extra \{dof\} the behavior is different: the output "my\_pressure" as defined previously is still correct around 2 having set also 2 in the fix command for the pressure\. But the temperature \(output as temp\) is now lower than the value set in fix, it is about 1\.68 where the value in the input is set to 2\. Other properties such as enthalpy are different\. By the way, the nr\. of dof to subtract is now "variable dof equal count\(spheroid\)" so one for each ellipsoid\. However I am puzzled by which is the \_real\_ temperature and pressure of the system, I mean the temperature and pressure that I should compare with a, say, NPT MonteCarlo simulation of the same GayBerne model\. Is the real pressure of my simulation that one in the default variable "press" or that one defined in the new variable "my\_pressure"? And For the temperature? Now is the default variable "temp" that appears scaled with respect to the value set in the input fix\.
If it is of any help I add here the input and a few lines from the output file\.
Thank you,
Giacomo
\#INPUT
units lj
atom\_style ellipsoid
dimension 3
lattice sc 0\.015
region box block 0 13 0 13 0 13
create\_box 1 box
create\_atoms 1 box
set type 1 mass 1\.0
set type 1 shape 4\.4 1 1
set group all quat 1 0 0 0
velocity all create 2\.4 87287 loop geom
pair\_style gayberne 1\.0 1\.0 1\.0 5\.5
pair\_coeff 1 1 1\.0 1\.0 2\.31 2\.31 0\.1155 0 0 0
neighbor 0\.8 bin
compute rot all temp/asphere
group spheroid type 1
variable dof equal count\(spheroid\)
compute\_modify rot extra {dof}
compute my_pressure all pressure rot
thermo_style custom step temp c_rot epair enthalpy etotal c_my_pressure vol
thermo 100
timestep 0.0015
compute q all property/atom quatw quati quatj quatk
dump 1 all custom 1000 trj.GB &
id type x y z c_q[1] c_q[2] c_q[3] c_q[4]
fix 1 all npt/asphere temp 2.00 2.00 0.1 iso 2.0 2.0 1
compute_modify 1_temp extra ${dof}
run 1000000
#OUPUT
#Step Temp rot E_pair Enthalpy TotEng my_press Volume
0 2.4 1.4393446 1763.3429 47315.517 1766.9413 683.22864 146466.67
1000 1.8456877 1.9774037 -0.82787275 5.2826365 1.9393986 0.24717927 42611.003
2000 1.8321838 1.9924696 -1.9225047 7.0596984 0.82452006 0.78871167 21517.332
3000 1.7931872 2.0044308 -2.8176202 9.9053969 -0.12906363 1.6336256 15577.271
4000 1.7155717 1.9612068 -3.1174212 11.384342 -0.54523505 2.0820965 14227.569
5000 1.6680423 2.0031676 -3.358239 10.57382 -0.85731451 2.0697931 13907.726
6000 1.6492788 1.9985652 -3.5217636 10.153491 -1.0489714 2.0529785 13788.782
7000 1.6903912 2.0188311 -3.622997 9.2808748 -1.0885643 1.9390379 13646.895
8000 1.68612 2.027784 -3.7486564 8.7942613 -1.2206277 1.8995213 13542.918
9000 1.6984239 2.0211098 -3.7926094 10.006982 -1.2461332 2.1108169 13451.657
10000 1.683795 1.9848341 -3.9342241 8.4912183 -1.4096812 1.896334 13353.377
11000 1.6507503 1.9869142 -4.0454844 8.7944424 -1.5704861 1.997618 13226.812
12000 1.6655222 2.0162935 -4.0906439 9.1784912 -1.5934978 2.0783132 13179.711
13000 1.656945 2.0052626 -4.1888793 8.8078819 -1.7045931 2.0447026 13106.972
14000 1.7271592 2.036114 -4.3834554 8.4367854 -1.7938958 2.0159741 12966.206
15000 1.6311151 1.9634758 -4.382641 8.603986 -1.9370819 2.0768644 12887.886
16000 1.6516237 1.9917753 -4.5394317 8.2047957 -2.0631239 2.033835 12894.296
17000 1.6451489 1.9738279 -4.7063156 8.2079605 -2.2397156 2.086091 12735.924
18000 1.6935545 2.0027681 -4.6540742 8.1704667 -2.1148987 2.0581861 12735.181
19000 1.6877883 1.9950598 -4.8879065 6.9418537 -2.3573763 1.9062623 12605.437
20000 1.696173 1.9953703 -4.9365715 7.7426354 -2.39347 2.06978 12489.544
21000 1.6338374 1.97224 -4.9679475 7.226103 -2.5183069 2.0010555 12514.519
22000 1.682566 2.001123 -5.0405691 6.4929179 -2.5178689 1.873131 12508.057
23000 1.6627053 2.0019965 -5.2120579 5.3724007 -2.7191351 1.7370322 12352.36
24000 1.7046418 2.0263275 -5.2014398 7.0500639 -2.645641 2.0331908 12285.028
25000 1.6375611 1.9730245 -5.3228392 7.0919864 -2.8676155 2.0925939 12190.478
26000 1.6563716 1.9732009 -5.4686041 5.6245524 -2.9851775 1.8448889 12197.693
27000 1.6714071 1.9564441 -5.5498537 7.4650911 -3.0438842 2.2214394 11965.96
28000 1.6674009 2.0019891 -5.6931811 6.338747 -3.1932182 2.0607704 11942.491
29000 1.7145486 2.0188442 -5.6585196 6.0503848 -3.0878673 1.9894478 11914.827
30000 1.6914707 2.0309485 -5.8145963 6.0237381 -3.2785452 2.0441763 11818.6
31000 1.652987 1.9988056 -5.9385707 6.2052696 -3.4602188 2.1310084 11695.908
32000 1.6480217 1.9813369 -6.0466228 5.3848725 -3.5757155 1.9896004 11722.128
33000 1.6391665 1.9509635 -5.7966 5.7619843 -3.3389693 1.9744756 11921.628
34000 1.7186204 2.0128303 -6.1053421 4.8597892 -3.5285849 1.8856257 11680.728
35000 1.6612858 1.9944786 -6.2583426 5.2401154 -3.7675481 2.0337424 11527.895
36000 1.7058402 2.0033237 -6.2068489 5.698711 -3.6492532 2.0933666 11525.433
37000 1.7288625 2.0616525 -6.1573918 5.433632 -3.5652785 2.0187548 11652.286
38000 1.6656828 1.9975877 -6.1869259 4.8016531 -3.689539 1.9235156 11599.495
39000 1.6507367 1.9985603 -6.2317657 5.5611394 -3.7567877 2.0988524 11513.21
40000 1.6791686 1.9943822 -6.2131265 6.151739 -3.69552 2.1976729 11489.306
41000 1.659162 2.0041769 -6.3180501 4.4441115 -3.8304399 1.9001133 11512.102
42000 1.6992902 2.0305375 -6.3068939 5.1836438 -3.7591188 2.0208522 11554.907
43000 1.6762376 1.9991044 -6.3031476 6.4593257 -3.7899356 2.2832563 11455.859
44000 1.713829 2.0449371 -6.2835754 6.0210137 -3.714002 2.1873334 11480.709
45000 1.7302254 2.0609315 -6.3160569 5.140519 -3.7219002 2.0222879 11480.86
46000 1.6418578 1.9805512 -6.3652749 4.5503099 -3.9036092 1.9347325 11484.719
I am puzzled by which is the _real_ temperature and pressure of the system,
I mean the temperature and pressure that I should compare with a, say, NPT
MonteCarlo simulation of the same GayBerne model. Is the real pressure of my
simulation that one in the default variable "press" or that one defined in
the new variable "my_pressure"? And For the temperature? Now is the default
variable "temp" that appears scaled with respect to the value set in the
input fix.
Here is the way to understand it. Fix npt/asphere will do everything
correctly (temp and pressure), except that the compute temp/asphere it
uses isn't smart enough to know about the uniaxial case Mike discussed. It
uses 6 dof per particle, instead of 5, as Mike says is needed for
uniaxial. Hence
you need to use compute modify extra N to remove one dof per particle.
And you must do this for the compute temp/asphere that fix npt/asphere is
using, as Mike indicated. Or any compute temp/asphere you create, e.g.
to monitor the temp.
Now if you want to monitor the temp and pressure that fix npt/asphere is
using, you can't just use the default "temp" and "press" in your thermo output.
As I said earlier, those just compute a temp for the 3 translational dof per
particle. You can do 1 of 2 things. You can assign the same computes that
fix npt/asphere is using to "temp" and "press" in thermo output, using
the thermo_modify
temp/press options. Or you can use thermo_style custom and access those
computes and output their values directly. Mike told you the names
of those computes, which are also listed on the fix npt doc page.
If you do that you should see an equilibrated temp and pressure that
is precisely what
you specified as targets for fix npt/asphere. Note however that if you
don't use compute modify as indicated above, then fix npt/asphere will
be computing the (slightly wrong) temp and thus pressure. Either way
you will still see the correct target values in your output, but if you make the
dof correction, the box volume should change (slightly) due to changing
the way the pressure will then be computed. I.e. a different box
size will give the same pressure as before.
Mike - I will add some comments to the compute temp/asphere doc
page about the uniaxial/biaxial issues.
Steve
How about this, to only be concerned with temperature including angular terms and pressure calculated using this temperature (1_temp is temperature, and 1_press is pressure):
#INPUT
units lj
atom_style ellipsoid
dimension 3
lattice sc 0.015
region box block 0 13 0 13 0 13
create_box 1 box
create_atoms 1 box
set type 1 mass 1.0
set type 1 shape 4.4 1 1
set group all quat 1 0 0 0
velocity all create 2.4 87287 loop geom
pair_style gayberne 1.0 1.0 1.0 5.5
pair_coeff 1 1 1.0 1.0 2.31 2.31 0.1155 0 0 0
neighbor 0.8 bin
thermo 100
timestep 0.0015
compute q all property/atom quatw quati quatj quatk
dump 1 all custom 1000 trj.GB &
id type x y z c_q[1] c_q[2] c_q[3] c_q[4]
fix 1 all npt/asphere temp 2.00 2.00 0.1 iso 2.0 2.0 1
# Count number of uniaxial particles
group uniaxial type 1
variable dof equal count(uniaxial)
# Remove a degree of freedom for each uniaxial particle from
# the NPT temperature compute
compute_modify 1_temp extra ${dof}
# Only output the temperature and pressure used by NPT/asphere
thermo_style custom step c_1_temp epair enthalpy &
etotal c_1_press vol
run 1000000
- Mike
Dear Steve and Mike,
yes, it looks fine now, thanks a lot for the detailed explanations and the input example.
Best wishes,
Giacomo