[gmx-users] Computing the energy of a crystal unitcell

[David van der Spoel]

On 2013-12-10 02:40, Golshan Hejazi wrote:
> Thanks for your reply.
>
> I have a paracetamol crystal and I would like to estimate the paracetamol/paracetamol interaction per unitcell.
>
> To do this I am doing the following:
>
> Energy of interaction of paracetamols per unitcell=
> Energy of a unitcell - (Energy of a molecule in the unitcell* number of the molecules in the unitcell)
>
> Energy of interaction of paracetamols per unitcell=  Term1 + Term2
>
>
> How to compute Term1:
> Took a crystal slab. Compute the total energy and divide it by the number of unitcells. The results is: -2443.29 kj/mol
>
> How to compute Term2:
> number of the molecules in the unitcell is 4. Energy of a single molecule in vacuum is 490.20 kj/mol
>
> Therefore:
> Energy of interaction of paracetamols per unitcell= -2443.29 kj/mol-(-490.20*4) kj/mol= -482.378 kj/mol
>
>
>
> Which is more than a covalent bond!
> Could you please help me to figure out what is going wrong?
>
The energy is only relative to something.

What you are computing (but the minus signs are a bit messy) is the 
interaction energy between the four monomers in a unit cell, which you 
still should divide by 4. In principle that should be comparable to the 
sublimation energy for paracetamol, if that were known, but -120 kJ/mol 
sounds reasonable. Why don't you try it for a simpler compound first?
> Thanks
> G.
>
>
>
>
> On Sunday, December 8, 2013 3:49 PM, David van der Spoel <spoel at xray.bmc.uu.se> wrote:
>
> On 2013-12-08 21:12, Golshan Hejazi wrote:
>> Hello everyone,
>>
>> I am trying to compute the energy of a crystal slab using steepest decent. Learning from previous questions in the group. I am using pbc in all direction. I need to compute this energy for one layer of crystal. Since crystal unit cell vectors (a b and c) can be less than rcoulomb and rlist ... it is obvious that I receive error while generating the tpr file.
>>
>> as an example, lets assume that a=0.7, b=.9 and c=1.1 angstrom, this means if I have a crystal slab of face(100), and b and c are in XY plane, then box dimension along z axes which is equal to 0.7 ... is less than rlist ...
>>
>> how can i compute the energy of a crystal slab then? If I decrease rcoulomb and rlist ... I am introducing artifact in computing electrostatic interaction and obviously, I can not increase box dimension along Z because in that way, crystal wont be periodic along Z anymore.
>>
>> Best and thanks for you useful comments in advance
>> G.
>>
> Use multiple layers and divide by the number of layers.
>


-- 
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se