towhee-users Mailing List for MCCCS Towhee

Pascal,
  Besides the options Ed mentioned, a fourth one would be to do this adsorption simulation in the
NPT ensemble using two boxes. The box where silicalite is present is kept fixed, on the other hand, 
pressure of the vapor box is regulated by performing volume exchange and particle swap moves. This
option will certainly relieve  you from all chemical potential related troubles, however, it is
computationaly a very expensive method. Catlow_Zeolite_4a example that comes in towhee package can
be a guide for you to set up such a simulation.

Ozgur

A. Ozgur Yazaydin
Post Doctoral Research Fellow
Department of Chemical & Biological Engineering
Northwestern University
2145 Sheridan Rd.
Evanston, IL 60208


Pascal Boulet wrote:
> Dear towhee users,
>
> I am currently facing a problem that puzzles me for a while.
>
> I have run GC simulations on the adsorption of para-cresol into 
> silicalite. The molecule size matches more or less exactely the size of 
> the zeolite channel. For getting at least one molecule adsorbed into the 
> zeolite I have to set the chemical potential to at least 200 kJ/mol 
> (note that the standard chemical potential equals to -75 kJ/mol). 
> Assuming that paracresol behaves like an ideal gas, the corresponding 
> pressure equals to about 1.0e17 Pa (if I am correct), which seems 
> meaningless (if not ridiculous!).
>
> I have seen in Snurr's and coll.'s paper (1993) that they simulated the 
> isotherm of adsorption of xylene in silicalite, a similar molecule to 
> paracresol. The equilibrium pressures reported are at most 1 kPa.
>
> How can I justify that the chemical potential must be that high in my 
> simulations? Is there anything wrong? Have you ever experienced such a 
> behaviour in your simulations?
>
> For information, I am using the Coupled to pre-nonbond algorithm and the 
> temperature is 310 K. The adsorption energies are in agreement with 
> experimental data.
>
> Thank you for your help
> Best Regards
> Pascal
>
>
>   

-------------------------------------------------------------------------
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse0120000070mrt/direct/01/_______________________________________________
Towhee-users mailing list
Tow...@li...
https://lists.sourceforge.net/lists/listinfo/towhee-users
===========End of original message text===========

Dear Randall and Ed, Dear all,

Thank you for your prompt reply. I come to the conclusion that I have a 
problem with the algorithm, in the sense that I am probably not using 
the right biasing (or the right parameters), as the acceptance ratio is 
very very small. Here's an example of towhee output that I get:

* Configurational-Bias REGROWTH Moves *
 Molecule type:    2 Box:    1
     Length  Attempts  Regrown  Accepted  %Regrown   %Accep.
        1      1571      1571       746    100.00     47.49
        3       761       761       669    100.00     87.91
       14      3890      3008         6     77.33      0.15
       15     93724     79086        32     84.38      0.03

The consequence is that, because of the too high chemical potential 
applied in the simulation, the calculated isotherm is probably *NOT* at 
the equilibrium. Indeed, there is nearly no fluctuation of the number of 
molecules during the course of the simulation (i.e., <N>= integer 
numbers at the end of the simulations).

To finish, Towhee is a great package: I did some nice simulations on the 
adsorption of CH4 into amorphous, porous silica, which admittedly is not 
a patological system. But Towhee has many interesting biased algorithms. 
And I am expecting to get good results with p-cresol by putting more 
efforts.

thank you again.
Best Regards
Pascal Boulet


-- 
Dr. Pascal Boulet, Computational Chemist
University of Provence
Laboratoire Chimie Provence, UMR6264
Centre Saint-Jerome, case MADIREL
F-13397 MARSEILLE Cedex 20, France
Tel. +33 (0) 491 63 71 17 Fax. +33 (0) 491 63 71 11
courriel: pas...@un...
http://www.lc-provence.fr
http://allos.up.univ-mrs.fr/boulet

Dear Pascal:

I assume you are interested in simulating the vapor adsorption of 
p-cresol. In this case, you must ensure that the pressure that 
corresponds to the chemical potential you have set does not exceed the 
vapor pressure of p-cresol at the temperature you are simulating. 
Otherwise, you are "simulating" a high pressure liquid in contact with 
silicalite. This is a big problem with specifying the chemical potential 
in a GCMC simulation - linking it back to the pressure is non-obvious. 
There are many examples in the literature where people have simulated 
vapor isotherms at chemical potentials that correspond to pressures well 
above the vapor pressure. They get interesting (but experimentally 
unattainable) results.

I think you have two options. You can relate chemical potential to 
pressure by performing an NPT simulation and computing the corresponding 
chemical potential with Widom insertion or some other technique, but 
this is not easy for such a big molecule. Or, you can explicitly relate 
the chemical potential to the fugacity and correct for any non-ideality 
with a fugacity coefficient. If you are treating p-cresol as a rigid 
molecule, then the ideal gas contribution to the chemical potential is 
easy to compute. If not, then you would have to evaluate the 
configurational integral of the molecule from a simulation of single 
molecule in the gas phase. The third option is to re-write the GCMC 
acceptance rules so that fugacity (and not chemical potential) are set 
directly. Marcus and I have had discussions about this, but this option 
is not yet supported in Towhee! The following paper may be helpful:

M. D. Macedonia and E. J. Maginn, “A Biased Grand Canonical Monte Carlo 
Method for
Simulating Adsorption Using All-Atom and Branched United Atom Models”, 
Molecular Physics, 1999, 96, 1375-1390.

Good luck,
Ed

Edward J. Maginn
Department of Chemical and Biomolecular Engineering
University of Notre Dame
Notre Dame, IN 46556 USA
email: ed...@nd...
internet: http://cbe.nd.edu
ph: 574. 631.5687
fax: 574.631.8366 



Pascal Boulet wrote:
> Dear towhee users,
>
> I am currently facing a problem that puzzles me for a while.
>
> I have run GC simulations on the adsorption of para-cresol into 
> silicalite. The molecule size matches more or less exactely the size of 
> the zeolite channel. For getting at least one molecule adsorbed into the 
> zeolite I have to set the chemical potential to at least 200 kJ/mol 
> (note that the standard chemical potential equals to -75 kJ/mol). 
> Assuming that paracresol behaves like an ideal gas, the corresponding 
> pressure equals to about 1.0e17 Pa (if I am correct), which seems 
> meaningless (if not ridiculous!).
>
> I have seen in Snurr's and coll.'s paper (1993) that they simulated the 
> isotherm of adsorption of xylene in silicalite, a similar molecule to 
> paracresol. The equilibrium pressures reported are at most 1 kPa.
>
> How can I justify that the chemical potential must be that high in my 
> simulations? Is there anything wrong? Have you ever experienced such a 
> behaviour in your simulations?
>
> For information, I am using the Coupled to pre-nonbond algorithm and the 
> temperature is 310 K. The adsorption energies are in agreement with 
> experimental data.
>
> Thank you for your help
> Best Regards
> Pascal
>
>
>   

Dear towhee users,

I am currently facing a problem that puzzles me for a while.

I have run GC simulations on the adsorption of para-cresol into 
silicalite. The molecule size matches more or less exactely the size of 
the zeolite channel. For getting at least one molecule adsorbed into the 
zeolite I have to set the chemical potential to at least 200 kJ/mol 
(note that the standard chemical potential equals to -75 kJ/mol). 
Assuming that paracresol behaves like an ideal gas, the corresponding 
pressure equals to about 1.0e17 Pa (if I am correct), which seems 
meaningless (if not ridiculous!).

I have seen in Snurr's and coll.'s paper (1993) that they simulated the 
isotherm of adsorption of xylene in silicalite, a similar molecule to 
paracresol. The equilibrium pressures reported are at most 1 kPa.

How can I justify that the chemical potential must be that high in my 
simulations? Is there anything wrong? Have you ever experienced such a 
behaviour in your simulations?

For information, I am using the Coupled to pre-nonbond algorithm and the 
temperature is 310 K. The adsorption energies are in agreement with 
experimental data.

Thank you for your help
Best Regards
Pascal


-- 
Dr. Pascal Boulet, Computational Chemist
University of Provence
Laboratoire Chimie Provence, UMR6264
Centre Saint-Jerome, case MADIREL
F-13397 MARSEILLE Cedex 20, France
Tel. +33 (0) 491 63 71 17 Fax. +33 (0) 491 63 71 11
courriel: pas...@un...
http://www.lc-provence.fr
http://allos.up.univ-mrs.fr/boulet

2008/2/23 Marc Segovia <ma...@fq...>:

>   I want to  simulate supercritical CO2,   with EMP2 forcefield, 32
>  molecules, using as reference (J.Chem.Phys.,120,9695), NPT simulation,
>  but I experienced some problems.
>
>  1) From the initial coordinates taken from pdb file, using pdb2towhee
>  utility, in the towhee start a new configuration appears, probably
>  related by center mass routine, ctrmas.F, this new configuration miss
>  most of molecules.

That frequently happens when you are starting a simulation of
molecules generated via some other program.  Towhee like to make sure
that all atoms in each molecule are in the same reference box (if
possible), and that the center-of-mass of each molecule is in the
central simulation box.  Just book-keeping issues, but important later
on for some of the moves.

>  2) and probably related with part 1), a message  about vibration
>  overlap, but this condition appears only with EMP2, not with others
>  forcefields like OPLS.

Most likely because OPLS uses flexible bond lengths, while EMP2 uses
rigid bond lengths.

>  Bond Types
>  Type:      2 Style: Fixed  Length: 1.1490

>  Initial Energies for Box     1
>   overlap vibration
>   ichain,imolty,ii,ip1 1 1 2 1
>   MAINLOOP: overlap in initial configuration

Towhee is detecting a bond length that is out of the tolerance for
rigid bond lengths and that is a fatal error.

Your best bet is to just generate a new simulation box from scratch.
CO2 is not especially difficult to equilibrate.  Just use the number
of molecules and box dimensions from your previous file and start
anew.

Marcus

-- 
Marcus G. Martin
Director, Useful Bias Incorporated
88 Martinez Road
Edgewood NM 87015-8222
(505) 286-4457
www.usefulbias.com
www.photobirder.com

Hello

 I want to  simulate supercritical CO2,   with EMP2 forcefield, 32 
molecules, using as reference (J.Chem.Phys.,120,9695), NPT simulation,  
but I experienced some problems.

1) From the initial coordinates taken from pdb file, using pdb2towhee 
utility, in the towhee start a new configuration appears, probably 
related by center mass routine, ctrmas.F, this new configuration miss 
most of molecules.

2) and probably related with part 1), a message  about vibration 
overlap, but this condition appears only with EMP2, not with others 
forcefields like OPLS.

Attached are pdb, towhee_coords, input and output file.

I would appreciate your help.

Best Regards and Nice weekend


Marc


Marc Segovia,PhD
Assistant Professor
Universidad de la Republica
Montevideo Uruguay
www.ccbg.fq.edu.uy

Hello,

This is far from my area, but a friend of mine who has done some work on 
adsorption has pointed to this paper:

Belfort, G. Similarity of Ideal Adsorbed Solution and
Potential theories for adsorption from a bulk phase onto a solid
surface. AIChE J. 1981, 27, 1021

Hopefully they'll discuss the details of IAST there.

Best,

Daniel


> What is the correct assumption? Do you know of any references?
>
> Best regards,
> Niels



-- 
Daniel Duque
http://rincon.uam.es/dir?cw=950067138671875

Dear towhee users,

a little bit off topic but still:

What are the underlying assumptions (on a molecular level) of the  
Ideal Adsorption Solution (IAS) theory introduced by Myers and  
Prausnitz in 1965?

While there is usually no disagreement on the following assumptions,
- equal accessibility of all adsorption sites to all adsorbates
- reversibility of adsorption
- perfect mixing

I found two different statements about adsorbate-adsorbate interactions:
(1) They are assumed to be equal for all adsorbates (A<->A = A<->B = B<->B)
(2) They are assumed to be equal to zero.

What is the correct assumption? Do you know of any references?

Best regards,
Niels

Dear Sir/Madam,

I am a PhD student of Hamburg University of Technology and working on
adsorption on zeolites. At this moment I am working on treating some
experimental data for NaX structure zeolite (with 86 Na atoms). I  
tried to simulate some data taken from S. Calero et al. .J. Am. Soc.  
2004, 126, p.113-11386 using Towhee codes. I am interested in the  
possibility for blocking
the sodalite cages within NaX structure. I have introduced dummy  
atoms, and made several trials with different sigma values (keeping an  
epsylon value = 0), and different values for rcutin. But the results  
of my simulations (on ethane adsorption) show very high values for  
ethane loading. It seems that rcutin does not have an impact when  
introducing dummy atoms. In attachment you can find input, forcefield  
and result files - 6 files: 3 force field files, 1 file with some more  
simulation results and 2 input files.

Could you please to give me your comments on it and describe the
possibility for blocking of sodalite cages using Towhee codes.

I am looking forward for your reply.

With best regards,

Diana Khashimova

Dpil.-Ing.

    Hamburg University of Technology
    Institute of Chemical Reaction Engineering
    Eisendorfer Str. 38
    21073 Hamburg, Germany
    Tel: +49 (0) 40 42878 2544
    Fax: +49 (0) 40 42878 2145

On Feb 2, 2008 1:52 PM, Craig Tenney <ct...@um...> wrote:

> I believe that for a non-rectangular simulation cell with periodic boundaries, the minimum image routine used in the twh_mimage function in mimage.f is guaranteed to return the true minimum image distance (rijtrue) only if rijtrue is less than one-half the minimum box width (Lmin/2) (ref 1). This can lead to errors when calculating two-body energies, as detailed below.

That is quite possible as I believe Towhee has infrequently been used
for anything other than cubic boxes.  Most of that use comes from
almost rectangular zeolite structures.

It would be very helpful if you could submit some input files that
(relatively quickly) generate the error.  Then I can work on debugging
the problem.

The best place to submit those files is the towhee-bugs mailing list.

thanks,

Marcus

-- 
Marcus G. Martin
Director, Useful Bias Incorporated
88 Martinez Road
Edgewood NM 87015-8222
(505) 286-4457
www.usefulbias.com
www.photobirder.com

Hello,

I believe that for a non-rectangular simulation cell with periodic boundaries, the minimum image routine used in the twh_mimage function in mimage.f is guaranteed to return the true minimum image distance (rijtrue) only if rijtrue is less than one-half the minimum box width (Lmin/2) (ref 1). This can lead to errors when calculating two-body energies, as detailed below.

twh_eng_twobody (in energy.f) and twh_engtotal_twobody (in engtotal.f) first test whether the center-of-mass (com) separation distance between two molecules exceeds the cutoff value rcm (rcm = rcelect + rcmi + rcmj for coulombic systems or rcm = rcut + rcmi + rcmj for vdw systems, where rcmx is the distance from the com of molecule x to its most distant atom). The test "if ( rijsq .gt. rcmsq ) lcompute = .false." will incorrectly ignore some two-body interactions between molecules when rcm exceeds Lmin/2 and twh_mimage unluckily returns an incorrect result.

This issue was first noticed because running and total energies did not agree at the end of some of my simulations. The disagreement was more significant for coulombic energies and simulations with larger molecules.

Regarding fixes I've tried or thought of: If twh_engtotal_twobody only runs at the start and finish of the simulation, turning off the "( rijsq .gt. rcmsq )" test provides a reliable check at the end of the simulation with minimal added computation. If there are problems, turning off the test in twh_eng_twobody will make them go away, generally at the cost of greater computing time. If possible, choosing a 'more rectangular' simulation cell sometimes helps. A more robust (and expensive) minimum image algorithm is possible (ref 2), but I've not tried implementing it. I've also not yet tried using 'ewald_fixed_cutoff' and setting rcelect and rcut less than (Lmin/2 - rcmi -rcmj), but I suspect it will work. It occurs to me now that dumping Lmin/2 and rcmx values to Towhee's output at the start of the simulation would make checking for this potential problem a relatively easy proposition.

Anyway, I'd not seen mention of this issue here before, so I figured documenting it might save other folks some debugging time in the future.

Regards,
Craig

References:
1. W. Smith, "The Minimum Image Convention in Non-Cubic MD Cells", Information Quarterly for Computer Simulation of Condensed Matter, No 30, April 1989.
2. M. Mezei, "Determining nearest image in non-orthogonal periodic systems", Information Quarterly for Computer Simulation of Condensed Phases, No 34, 1992.