[gmx-users] Three-site model for acetonitrile

[Christoph Freudenberger]

Anton Feenstra wrote:
> It seems quite OK, however it seems you simply moved mass from the
> C to the Me and N atoms (I'm not sure, since I don't know the exact
> masses of the atoms). This will influence the moments of intertia
> of your molecule (i.e. it will tumble slower).
> 
> The way to do this more accurately, is to calculate the inertia tensor
> for your normal, three mass molecule, and then calculate the distance
> and masses for a two mass system with the same total mass and the same
> inertia tensors. For a linear molecule, this is quite trivial to do.
> You will find that the two masses will need to be closer together than
> your Me and N atoms are. Then, if you have the two masses at their
> appropriate distance, you construct *all* *three* atoms (Me, C and N)
> from the two dummy masses.
> 
> You will now have a system of two masses that accurately reproduce
> the mass and moments of inertia of acetonitile and three 'dummy' atoms
> that accurately reproduce the normal atomic interactions of acetonitrile.
Thanks for your suggestions... yes I just moved the mass of the carbon
into the methyl and N with the condition of the com to be preserved.

setting dispcorr improved things a little bit but the densities are
still too low:
expt   dispcorr=no  dispcorr=enerpres
780       720            740

So maybe it's due to the error in the moment of inertia...

Just to see if I did every thing right:
The com of MeCN I get from:
M = sum(mi)
rcom = sum(mi*ri)/M
The inertia tensor of a linear molecule is:
     [Ia 0 0]
I = [0 Ia 0] with Ia=sum(mi(ri-rcom)^2)
     [0  0 0]
Now I have to choose to masses m1 and m2 with r1 and r2 relativ to the com:
The distance R between m1 and m2 can be chosen freely:
R = r1 - r2
There are three conditions to be satisfied:
M = m1 + m2
m1*r1 + m2*r2 = 0
m1*r1^2 + m2*r2^2 = Ia
Using these 4 eqn I can calculate the four unknows.

The new topology would come out this way:
;------------------------------------------------------
[ atomtypes ]
;type   mass           charge    ptype c6            c12
D1       9.49031       0.000       D   0.0           0.0
D2      31.56239       0.000       D   0.0           0.0
MeAN     0.00000       0.000       A   0.90571E-02   0.26212E-04 ; (*)
  CAN     0.00000       0.000       A   0.51454E-02   0.12167E-04 ; (*)
  NAN     0.00000       0.000       A   0.26955E-02   0.28943E-05 ; (*)

[ moleculetype ]
; name  nrexcl
Acetonitril        2

[ atoms ]
;   nr  type    resnr   residu  atom    cgnr    charge  mass
1       D1       1       MeCN    D1AN    1        0.000   9.49031
2       D2       1       MeCN    D2AN    1        0.000  31.56239
3       MeAN     1       MeCN    C3AN    1        0.206   0.00000
4       CAN      1       MeCN    C4AN    1        0.247   0.00000
5       NAN      1       MeCN    N5AN    1       -0.453   0.00000

[ constraints ]
;  ai  aj funct           b0
1       2       1          0.26300

[ dummies2 ]
;  ai    aj    ak       funct   a
     3     1     2       1       0.2652197
     4     1     2       1       0.8203527
     5     1     2       1       1.2652197

[ exclusions ]
3       4       5
4       5       3
5       4       3

;(*) E.Guardia, Mol. Simul. 26, 2001, 2878;
;A rigid two mass system with the same mom. of inertia
;as MeCN is used to provide a linear molecule. Me, C and
;N are defined as dummies relativ to D1 and D2.
;---------------- eof -------------------

How does that look...? I will start the new simulation
as soon as possible

regards
-- 
Christoph Freudenberger
Abt. Organische Chemie I AK Siehl - Uni Ulm -Tel: ++49-731-502-2785