I am studying the formation of a carbon-based porous material using a coarse-grained ellipsoid model and the gayberne potential for anisotropic interactions.
As suggested from the manual/examples/mailing list I am using:
compute mytemp all temp/asphere
with fix rigid/nve/small + fix langevin. (for brownian dynamics)
This gives an error "Temperature compute degrees of freedom < 0" because the default temperature (translations only) compute is being used, hence I override this using:
thermo_modify temp mytemp.
My desired temperature for fix langevin is 1.0 in LJ units (NOTE: the temperature used by this fix is set to mytemp using fix_modify) but the thermo (set to c_mytemp) output is showing equilibration at a different temperature (1.5 in this case).
I am using the latest version of LAMMPS (7-Dec-15) and the angmom option on fix langevin is off - I have tried with angmom turned on and set to 1 or 3.33 as suggested but this also gives the incorrect temperature and in some cases it causes lost atoms.
Below is an excerpt from my input file that deals with the temperature settings:
What is the cause of this issue? I believe it is likely the amount of degrees of freedom being taken into account, but compute temp/asphere should automatically subtract all degrees of freedom until there are 6 per rigid body (correct?) so I remain quite confused.
It is important to note that I do not have this problem when I simulate with fix nve/asphere (no rigid bodies are present, hence no fix rigid) and set the angmom factor to 3.33 in fix langevin.
I have spent a long time on the mailing list and the manual searching for a solution and haven't found one but there is always a chance that I missed something hence I have posted here, but hopefully I am not wasting anyone's time.
I don’t understand the rigid body part of your model. If you
are modeling a collection of ellipsoids with the GB potential,
what are the rigid bodies which would require use of fix rigid/nve?
The ellipsoids themselves are not rigid bodies in LAMMPS lingo.
They are just individual aspherical particles.
There is a very specific case in which compute temp/asphere seems to be inappropriate for calculating the temperature of a system of ellipsoids undergoing langevin dynamics (using fix langevin), and that is for systems of rigid bodies made up of ellipsoidal particles. The cause of this I have not been able to determine, but it is important to note that the kinetic energy (from thermo output or compute ke + compute erotate/asphere) implies that the number of degrees of freedom are being calculated correctly (KE = (3N -3)kt , N = number of rigid bodies, 6 DOF per rigid body).
Using fix rigid/small langevin and the temperature output by this fix (rigid/nve/small doesn't output the temperature and it isn't highlighted in the manual exactly what output value corresponds to) everything works as it should. Note the DOF are the same and the comparison with KE (which is = compute erotate/rigid + compute ke/rigid) gives the same result as that mentioned above, ie. the DOF are calculated correctly as 6 per rigid body.
From looking at the code I have not been able to determine the cause of the issue and perhaps I am missing something but I will highlight that fix langevin coupled with rigid bodies of LJ spheres works (compute temp - default) and fix langevin coupled with free or bonded (but not rigid) ellipsoids (compute temp/asphere) also works.
Attached is an input script, a data file and a log file that will highlight the issue as the thermo output contains all values I mentioned above (Everything is in LJ units). I believe there should be a note in the manual highlighting this specific case for future users and I am more then happy to include more test cases or discuss further with the developers to help fix this issue.
Averages of all the columns in the example log file (Langevin desired temperature = 1.0):
rottemp = 1.52776899 ± 0.125528157205
KinEng = 0.907494773 ± 0.0745637155198 (Total = 226.87369325 [* no_atoms])
asphke = 0.907494773 ± 0.0745637155198
rigid = 1.023166451 ± 0.079194023097
rigtotke = 0.613899879 ± 0.0475164140354 (Total = 153.47496975 [* no_atoms])
DOF = 147
Finally, it is necessary for my future work that I can set multiple damping parameters for different types of atoms - is there a reason this hasn't been implemented in the fix rigid langevin code?
Regards,
Andrew
Andrew Tarzia. Doctoral Candidate.
School of Chemistry and Physics | University of Adelaide | [email protected]... |
Thank you Steve, you have interpreted everything I said correctly which is sometimes hard to do.
Is it true that the use of multiple fix rigid commands (to obtain different damping factors for different rigid bodies as you suggested) greatly effects performance? There is a warning on the doc page for fix rigid against doing this hence I was aiming to avoid it.
I should be more clear and say that in the future I would like to run simulations with multiple types of rigid bodies and hence multiple damping factors.
Thanks,
Andrew
Andrew Tarzia. Doctoral Candidate.
School of Chemistry and Physics | University of Adelaide | [email protected]... |
Is it true that the use of multiple fix rigid commands (to obtain different damping factors for different rigid bodies as you >suggested) greatly effects performance? There is a warning on the doc page for fix rigid against doing this hence I was aiming >to avoid it.
Hard to say. You’d have to time it for your problem.
I should be more clear and say that in the future I would like to run simulations with multiple types of rigid bodies and hence >multiple damping factors.
How are you proposing the fix rigid command would be able to distinguish
different kinds of rigid bodies? How would you label a rigid body
and tell the command to apply one damping factor to one label,
and another damping factor to another label? That is what
atom types do for fix langevin. There is no such “type” for a rigid body.
Perhaps knowing the types of atoms in each rigid body, assuming that type 1 and 2 for example were exclusive to one type of rigid body, and then applying fix rigid to a group that contains those atoms? That seems like the most logical way to do it off the top of my head...
Andrew
Andrew Tarzia. Doctoral Candidate.
School of Chemistry and Physics | University of Adelaide | [email protected]... |