GCMC
395 views
Skip to first unread message
Isaak Daniels
Aug 20, 2011, 1:51:27 PM8/20/11
to cp2k
Is there a way to run a sort of GCMC on Cp2k?
Thank you
Thank you
Matt McGrath
Aug 21, 2011, 8:42:42 PM8/21/11
to cp2k
Hi Isaak. Unfortunately, this isn't really possible at the moment,
unless you do Gibbs ensemble Monte Carlo and make the second box huge,
i.e. essentially a reservoir. Of course, if you're using any grid-
based method to compute the energy (Quickstep, Ewald) this adds a
whole lot of expense to the simulation as well, and probably isn't
worth it.
GCMC has been in the pipeline for a while (there have been issues with
the choice of reference state for QS calculations with GCMC), and is
the next big project on my list. I'm about 30% done with the current
project (making the MC routines more streamlined), so it might get
implemented in the next year...but that's just a guess. Could be six
months, could be a couple years, so if you're working on a major
project, it's probably best to use a different code. Sorry!
Cheers, Matt
unless you do Gibbs ensemble Monte Carlo and make the second box huge,
i.e. essentially a reservoir. Of course, if you're using any grid-
based method to compute the energy (Quickstep, Ewald) this adds a
whole lot of expense to the simulation as well, and probably isn't
worth it.
GCMC has been in the pipeline for a while (there have been issues with
the choice of reference state for QS calculations with GCMC), and is
the next big project on my list. I'm about 30% done with the current
project (making the MC routines more streamlined), so it might get
implemented in the next year...but that's just a guess. Could be six
months, could be a couple years, so if you're working on a major
project, it's probably best to use a different code. Sorry!
Cheers, Matt
Isaak Daniels
Aug 22, 2011, 9:57:39 AM8/22/11
to cp2k
Can one do it so that the box not of interest (representing the
environment) is MM while the other is done with QM, so one can avoid
the pitfalls of enlarging the box?
Thank you
environment) is MM while the other is done with QM, so one can avoid
the pitfalls of enlarging the box?
Thank you
Matt McGrath
Aug 23, 2011, 4:21:41 AM8/23/11
to cp2k
Hmm. I was thinking about that, and I didn't think it would work
before I responded to your original question. But now, I'm not so
sure. If you set all the MM interactions in the vapor box to have
epsilon=0 and no charges (so you don't need Ewald sums), then you have
an ideal gas, so...maybe.
A quick way to check would be to try it with MM in both boxes, but set
all the LJ parameters and charges equal to zero in the vapor box.
Choose it to be large enough that the density doesn't really change
over the course of the simulation, don't do any molecule displacements
in the vapor box, and compare it against GCMC with a different code.
One concern I have is that, on swap moves, the code is going to check
for overlaps in the vapor box, which it should do for true GCMC. If
you try this test and get an answer that is close, and you feel
comfortable changing the source code, I can let you know which lines
to change to remove those checks so we can see if it works. I know of
codes that do GCMC as basically GEMC with a box that has no
interactions in it (and is big...500 molecules or more), so this
doesn't seem so different.
Cheers, Matt
On Aug 22, 10:57 pm, Isaak Daniels <isaakdanielsfi...@gmail.com>
wrote:
before I responded to your original question. But now, I'm not so
sure. If you set all the MM interactions in the vapor box to have
epsilon=0 and no charges (so you don't need Ewald sums), then you have
an ideal gas, so...maybe.
A quick way to check would be to try it with MM in both boxes, but set
all the LJ parameters and charges equal to zero in the vapor box.
Choose it to be large enough that the density doesn't really change
over the course of the simulation, don't do any molecule displacements
in the vapor box, and compare it against GCMC with a different code.
One concern I have is that, on swap moves, the code is going to check
for overlaps in the vapor box, which it should do for true GCMC. If
you try this test and get an answer that is close, and you feel
comfortable changing the source code, I can let you know which lines
to change to remove those checks so we can see if it works. I know of
codes that do GCMC as basically GEMC with a box that has no
interactions in it (and is big...500 molecules or more), so this
doesn't seem so different.
Cheers, Matt
On Aug 22, 10:57 pm, Isaak Daniels <isaakdanielsfi...@gmail.com>
wrote:
Isaak Daniels
Aug 23, 2011, 8:23:17 AM8/23/11
to cp2k
Can I also have it so that in the non-vapor box, gas molecules can
enter that box and one can also import coordinates of a crystal, with
this box having QM?
Thank you
enter that box and one can also import coordinates of a crystal, with
this box having QM?
Thank you
Matt McGrath
Aug 23, 2011, 11:08:15 PM8/23/11
to cp2k
When you say "import coordinates of a crystal", do you mean as a swap
move? As in, you want to have gas molecules and crystal molecules
come into the box, or that you put crystal coordinates into the box at
the beginning, treat the whole box as QM, and then swap gas molecules
into the box (also treated as QM)? The latter option should be
possible if the test I described above worked (since the force_env of
each box can be specified separately), and is what I assumed you
wanted to do in the first post. If you also want to swap in crystal
nuclei, you could do that by making it a new molecule type (and
creating a psf file for it), but then it would swap into the QM box in
a completely random position, and would still be treated as a molecule
from CP2K's standpoint (rotated and translated together). I don't
know if that's what you want or not.
Cheers, Matt
move? As in, you want to have gas molecules and crystal molecules
come into the box, or that you put crystal coordinates into the box at
the beginning, treat the whole box as QM, and then swap gas molecules
into the box (also treated as QM)? The latter option should be
possible if the test I described above worked (since the force_env of
each box can be specified separately), and is what I assumed you
wanted to do in the first post. If you also want to swap in crystal
nuclei, you could do that by making it a new molecule type (and
creating a psf file for it), but then it would swap into the QM box in
a completely random position, and would still be treated as a molecule
from CP2K's standpoint (rotated and translated together). I don't
know if that's what you want or not.
Cheers, Matt
Isaak Daniels
Aug 24, 2011, 8:21:22 AM8/24/11
to cp2k
Is there a way to set a chemical potential, such as -454?
Thank you
Thank you
Matt McGrath
Aug 25, 2011, 3:23:36 AM8/25/11
to cp2k
The only way to set the chemical potential is if you use the above
tricks and they work. You can't set it directly, but you could choose
your ideal gas to have a certain chemical potential. GEMC would
equilibrate the chemical potential between the two boxes, and provided
your ideal vapor box is big enough (i.e. a true reservoir), its
chemical potential would not change, so therefore the chemical
potential in both boxes would be (approximately) the same as what you
started with for your ideal gas. Finding what is "big enough" would
probably take some effort, though, and there would be error bars
associated with the number, so this is not a perfect solution.
Cheers, Matt
tricks and they work. You can't set it directly, but you could choose
your ideal gas to have a certain chemical potential. GEMC would
equilibrate the chemical potential between the two boxes, and provided
your ideal vapor box is big enough (i.e. a true reservoir), its
chemical potential would not change, so therefore the chemical
potential in both boxes would be (approximately) the same as what you
started with for your ideal gas. Finding what is "big enough" would
probably take some effort, though, and there would be error bars
associated with the number, so this is not a perfect solution.
Cheers, Matt
Isaak Daniels
Aug 26, 2011, 12:12:50 PM8/26/11
to cp2k
Is there a way to find the average number of molecules in each box for
a GEMC simulation?
Thank you
a GEMC simulation?
Thank you
Isaak Daniels
Aug 26, 2011, 4:07:39 PM8/26/11
to cp2k
And one would set it by using the Pressure?
Thank you
On Aug 26, 11:12 am, Isaak Daniels <isaakdanielsfi...@gmail.com>
wrote:
Thank you
On Aug 26, 11:12 am, Isaak Daniels <isaakdanielsfi...@gmail.com>
wrote:
> Is there a way to find the average number of molecules in each box for
> a GEMC simulation?
>
> Thank you
>
> On Aug 25, 2:23 am, Matt McGrath <obfisca...@gmail.com> wrote:
>
>
>
> > The only way to set thechemical potentialis if you use the above
> a GEMC simulation?
>
> Thank you
>
> On Aug 25, 2:23 am, Matt McGrath <obfisca...@gmail.com> wrote:
>
>
>
> > tricks and they work. You can't set it directly, but you could choose
> > your ideal gas to have a certainchemical potential. GEMC would
> > equilibrate thechemical potentialbetween the two boxes, and provided
> > your ideal gas to have a certainchemical potential. GEMC would
> > your ideal vapor box is big enough (i.e. a true reservoir), its
> >chemical potentialwould not change, so therefore thechemical> potentialin both boxes would be (approximately) the same as what you
> > started with for your ideal gas. Finding what is "big enough" would
> > probably take some effort, though, and there would be error bars
> > associated with the number, so this is not a perfect solution.
>
> > Cheers, Matt
>
> > On Aug 24, 9:21 pm, Isaak Daniels <isaakdanielsfi...@gmail.com> wrote:
>
> > > Is there a way to set achemical potential, such as -454?
> > probably take some effort, though, and there would be error bars
> > associated with the number, so this is not a perfect solution.
>
> > Cheers, Matt
>
> > On Aug 24, 9:21 pm, Isaak Daniels <isaakdanielsfi...@gmail.com> wrote:
>
Isaak Daniels
Aug 26, 2011, 5:10:15 PM8/26/11
to cp2k
Also, how does one make sure the crystal is rigid but the atoms can go
in and out of each box?
Thank you for answering my questions
in and out of each box?
Thank you for answering my questions
Isaak Daniels
Aug 26, 2011, 5:26:59 PM8/26/11
to cp2k
Sorry for not being clear. By crystal, I mean a large molecule I want
rigid, and when I say atoms, I meant the gas atoms.
Thank you
rigid, and when I say atoms, I meant the gas atoms.
Thank you
Matt McGrath
Aug 28, 2011, 1:46:39 AM8/28/11
to cp2k
Hi Issak. You can find the average number of molecules in each box by
writing a script to average the mc_molecules file that is printed
out. The format is
Step number # of molecules of type 1 # of molecules of type
2 ..... etc, for box 1
Step number # of molecules of type 1 # of molecules of type
2 ..... etc, for box 2 (if there is a second box in the system)
At one time this was printed out at the end of a run, but for some
reason I removed it. There is no way to set this number via the
pressure...you just need to change the box volume, run the simulation,
see what the average density is (assuming the simulation has
equilibrated), change the box volume so that average density will give
the desired number of molecules, and then rerun. If the number of
molecules in the box is very small, there will be large fluctuations
around this number over the course of the run, but assuming the
simulation is long enough and equilibrated, the average should be what
you desire.
The easiest way of keeping a large molecule rigid is to change the
_MOL probabilities in the input file. For example, if the big
molecule is type 1 and the gas atoms are type 2, then a line of
PMSWAP_MOL 0.0 1.0
will let the gas molecules swap between boxes but not the big
molecule. Similarly,
PMTRAION_MOL 0.0 1.0
will attempt to change the conformation of molecule 2 but not molecule
1. If you want 30% of molecule translations to be done on molecule 1
and 70% on molecule 2, a line like
PMTRANS_MOL 0.3 1.0
would work. Hope that clarifies things. The file, tests/MC/regtest/
GEMC_NVT_box1.inp uses this, since it has two molecule types.
Cheers, Matt
writing a script to average the mc_molecules file that is printed
out. The format is
Step number # of molecules of type 1 # of molecules of type
2 ..... etc, for box 1
Step number # of molecules of type 1 # of molecules of type
2 ..... etc, for box 2 (if there is a second box in the system)
At one time this was printed out at the end of a run, but for some
reason I removed it. There is no way to set this number via the
pressure...you just need to change the box volume, run the simulation,
see what the average density is (assuming the simulation has
equilibrated), change the box volume so that average density will give
the desired number of molecules, and then rerun. If the number of
molecules in the box is very small, there will be large fluctuations
around this number over the course of the run, but assuming the
simulation is long enough and equilibrated, the average should be what
you desire.
The easiest way of keeping a large molecule rigid is to change the
_MOL probabilities in the input file. For example, if the big
molecule is type 1 and the gas atoms are type 2, then a line of
PMSWAP_MOL 0.0 1.0
will let the gas molecules swap between boxes but not the big
molecule. Similarly,
PMTRAION_MOL 0.0 1.0
will attempt to change the conformation of molecule 2 but not molecule
1. If you want 30% of molecule translations to be done on molecule 1
and 70% on molecule 2, a line like
PMTRANS_MOL 0.3 1.0
would work. Hope that clarifies things. The file, tests/MC/regtest/
GEMC_NVT_box1.inp uses this, since it has two molecule types.
Cheers, Matt
Isaak Daniels
Aug 28, 2011, 12:06:05 PM8/28/11
to cp2k
Just to confirm, there is no way to directly have an initial number of
molecules set that can be changed, for example doing 400 gas molecules
that can go in and out of boxes?
Thank you
molecules set that can be changed, for example doing 400 gas molecules
that can go in and out of boxes?
Thank you
Isaak Daniels
Aug 28, 2011, 1:36:32 PM8/28/11
to cp2k
Also, how do you do it so that one can specify a number of molecules
without needing coordinates?
Thank you for helping me.
On Aug 28, 11:06 am, Isaak Daniels <isaakdanielsfi...@gmail.com>
> > > > > > > > > > > in the vapor box, and compare it againstGCMCwith a different code.
> > > > > > > > > > > > > the choice of reference state for QS calculations withGCMC), and is
without needing coordinates?
Thank you for helping me.
On Aug 28, 11:06 am, Isaak Daniels <isaakdanielsfi...@gmail.com>
> > > > > > > > > > > One concern I have is that, on swap moves, the code is going to check
> > > > > > > > > > > for overlaps in the vapor box, which it should do for trueGCMC. If
> > > > > > > > > > > you try this test and get an answer that is close, and you feel
> > > > > > > > > > > comfortable changing the source code, I can let you know which lines
> > > > > > > > > > > to change to remove those checks so we can see if it works. I know of
> > > > > > > > > > > codes that doGCMCas basically GEMC with a box that has no
> > > > > > > > > > > comfortable changing the source code, I can let you know which lines
> > > > > > > > > > > to change to remove those checks so we can see if it works. I know of
> > > > > > > > > > > interactions in it (and is big...500 molecules or more), so this
> > > > > > > > > > > doesn't seem so different.
>
> > > > > > > > > > > Cheers, Matt
>
> > > > > > > > > > > On Aug 22, 10:57 pm, Isaak Daniels <isaakdanielsfi...@gmail.com>
> > > > > > > > > > > wrote:
>
> > > > > > > > > > > > Can one do it so that the box not of interest (representing the
> > > > > > > > > > > > environment) is MM while the other is done with QM, so one can avoid
> > > > > > > > > > > > the pitfalls of enlarging the box?
>
> > > > > > > > > > > > Thank you
>
> > > > > > > > > > > > On Aug 21, 7:42 pm, Matt McGrath <obfisca...@gmail.com> wrote:
>
> > > > > > > > > > > > > Hi Isaak. Unfortunately, this isn't really possible at the moment,
> > > > > > > > > > > > > unless you do Gibbs ensemble Monte Carlo and make the second box huge,
> > > > > > > > > > > > > i.e. essentially a reservoir. Of course, if you're using any grid-
> > > > > > > > > > > > > based method to compute the energy (Quickstep, Ewald) this adds a
> > > > > > > > > > > > > whole lot of expense to the simulation as well, and probably isn't
> > > > > > > > > > > > > worth it.
>
> > > > > > > > > > > > >GCMChas been in the pipeline for a while (there have been issues with
> > > > > > > > > > > doesn't seem so different.
>
> > > > > > > > > > > Cheers, Matt
>
> > > > > > > > > > > On Aug 22, 10:57 pm, Isaak Daniels <isaakdanielsfi...@gmail.com>
> > > > > > > > > > > wrote:
>
> > > > > > > > > > > > Can one do it so that the box not of interest (representing the
> > > > > > > > > > > > environment) is MM while the other is done with QM, so one can avoid
> > > > > > > > > > > > the pitfalls of enlarging the box?
>
> > > > > > > > > > > > Thank you
>
> > > > > > > > > > > > On Aug 21, 7:42 pm, Matt McGrath <obfisca...@gmail.com> wrote:
>
> > > > > > > > > > > > > Hi Isaak. Unfortunately, this isn't really possible at the moment,
> > > > > > > > > > > > > unless you do Gibbs ensemble Monte Carlo and make the second box huge,
> > > > > > > > > > > > > i.e. essentially a reservoir. Of course, if you're using any grid-
> > > > > > > > > > > > > based method to compute the energy (Quickstep, Ewald) this adds a
> > > > > > > > > > > > > whole lot of expense to the simulation as well, and probably isn't
> > > > > > > > > > > > > worth it.
>
> > > > > > > > > > > > > the choice of reference state for QS calculations withGCMC), and is
Matt McGrath
Aug 29, 2011, 1:44:59 AM8/29/11
to cp2k
Hi Issak. Neither of those is possible at the moment. The only
option is to write a script to generate the coordinates yourself (for
example, creating a grid inside a box the size of your desired
simulation box, and putting a molecule on each gridpoint, and then
using those coordinates as the input to one of the boxes in CP2K). If
you are using an ideal gas box (i.e. all the FIST interactions are set
to zero), you don't have to equilibrate the coordinates at all, which
is a small help. Otherwise you should run some short equilibration
(not letting the molecules swap) before doing a real run.
Cheers, Matt
option is to write a script to generate the coordinates yourself (for
example, creating a grid inside a box the size of your desired
simulation box, and putting a molecule on each gridpoint, and then
using those coordinates as the input to one of the boxes in CP2K). If
you are using an ideal gas box (i.e. all the FIST interactions are set
to zero), you don't have to equilibrate the coordinates at all, which
is a small help. Otherwise you should run some short equilibration
(not letting the molecules swap) before doing a real run.
Cheers, Matt
Isaak Daniels
Aug 29, 2011, 11:51:26 AM8/29/11
to cp2k
Thank you. Is there a way to stop the program from aborting when there
are no molecules left in one box?
Thank you
> ...
>
> read more »
are no molecules left in one box?
Thank you
>
> read more »
Matt McGrath
Aug 30, 2011, 3:41:17 AM8/30/11
to cp2k
Nope, not possible. I never saw a reason for having such a feature
(zero particles in the Gibbs ensemble vapor box is still a valid
configuration).
Cheers, Matt
On Aug 30, 12:51 am, Isaak Daniels <isaakdanielsfi...@gmail.com>
> > > > > > > > > > > > > > > > unless you do Gibbs ensemble Monte Carlo and make the second box huge,...
>
> read more »
(zero particles in the Gibbs ensemble vapor box is still a valid
configuration).
Cheers, Matt
On Aug 30, 12:51 am, Isaak Daniels <isaakdanielsfi...@gmail.com>
>
> read more »
Isaak Daniels
Aug 30, 2011, 11:14:35 AM8/30/11
to cp2k
What variables should one vary to affect the chemical potential?
Thank you
> ...
>
> read more »
Thank you
>
> read more »
Matt McGrath
Aug 31, 2011, 12:21:10 AM8/31/11
to cp2k
The chemical potential of an ideal gas is a function of the
temperature and number density (N/V). So, if you have a chemical
potential you want at a certain temperature, you can solve that
equation for the number density. Once you have the number density,
you decide how many vapor molecules you want (400, 500 maybe) and
calculate the vapor box volume from that. Then you run the simulation
and make sure the vapor density does not change over the course of the
simulation, because if it does, your chemical potential changes as
well.
Cheers, Matt
On Aug 31, 12:14 am, Isaak Daniels <isaakdanielsfi...@gmail.com>
> > > > > > > > > > > > > > > > interactions in it (and is...
>
> read more »
temperature and number density (N/V). So, if you have a chemical
potential you want at a certain temperature, you can solve that
equation for the number density. Once you have the number density,
you decide how many vapor molecules you want (400, 500 maybe) and
calculate the vapor box volume from that. Then you run the simulation
and make sure the vapor density does not change over the course of the
simulation, because if it does, your chemical potential changes as
well.
Cheers, Matt
On Aug 31, 12:14 am, Isaak Daniels <isaakdanielsfi...@gmail.com>
>
> read more »
Isaak Daniels
Sep 2, 2011, 3:24:16 PM9/2/11
to cp2k
How does one make sure the density isn't changed in the vapor box when
one moves a molecule out of it?
Thank you
> ...
>
> read more »
one moves a molecule out of it?
Thank you
>
> read more »
Matt McGrath
Sep 3, 2011, 1:56:43 AM9/3/11
to cp2k
Hi.
That's not possible. The more molecules you have, the less the
density will change, but it will always change by a little.
Cheers, Matt
> > > > > > > > > > > > > > > > > > before I responded to your original question. But now, I'm not so...
>
> read more »
That's not possible. The more molecules you have, the less the
density will change, but it will always change by a little.
Cheers, Matt
>
> read more »
Isaak Daniels
Sep 30, 2011, 2:02:12 PM9/30/11
to cp2k
I think I will try to implement GCMC into the source code. What files
should one modify if one wants to introduce a new ensamble?
Thank you
> ...
>
> read more »
should one modify if one wants to introduce a new ensamble?
Thank you
>
> read more »
Matt McGrath
Oct 1, 2011, 12:57:14 AM10/1/11
to cp2k
CP2K has a pretty consistent naming scheme. Usually, routines related
to a particular topic begin with TOPIC_. In this case, the Monte
Carlo ensemble routines are in mc_ensembles.F. From there you should
be able to figure out how the move types are selected (the definition
of an ensemble is just the selection of different move types).
Cheers, Matt
> > > > > > > > > > > > > > > > > > a completely random position, and would still be...
>
> read more »
to a particular topic begin with TOPIC_. In this case, the Monte
Carlo ensemble routines are in mc_ensembles.F. From there you should
be able to figure out how the move types are selected (the definition
of an ensemble is just the selection of different move types).
Cheers, Matt
>
> read more »
Isaak Daniels
Oct 1, 2011, 1:21:29 AM10/1/11
to cp2k
Thank you
Wish me luck
> ...
>
> read more »
Wish me luck
>
> read more »
huy pham
Jul 5, 2020, 12:11:17 PM7/5/20
to cp2k
Hi Isaak,
Best,
Did you have any luck implementing GCMC in CP2K?
I also want to run GCMC using CP2K for the MOF. Could you please share the code with me?
Thank you!
Best,
Huy
Reply all
Reply to author
Forward
0 new messages