[ase-users] Fw: MD simulation with NPT

[Gaël Donval]

Hi Shi,

Keep it into the list please by using "Reply All" instead of simply "Reply" (those lists are archived so that people can find issues and fixes in the future).

Yes, I was wrong about just swapping rows... You have to swap all rows and all columns so that:

C = [
    [a, 0, 0],
    [b, d, 0],
    [c, e, f],
]

becomes:

C' = [
    [f, e, c],
    [0, d, b],
    [0, 0, a],
]

(that is still in addition to transforming positions [x, y, z] to [z, y, x]).

Please note it is *not* a transpose but the application of a basis change to the matrix C...



Basically, we have that relationship between scaled (sx, ...) and cartesian coordinates (x, ...) where '@' is the dot product:

[sx, sy, sz] @ C = [x, y, z]

Let P be the matrix transforming any vector [u, v, w] into [w, v, u]:

[u, v, w] @ P = [w, v, u]

Then, using the identity matrix I = P @ P^{-1}:

[sx, sy, sz] @ I @ C = [x, y, z]
[sx, sy, sz] @ (P @ P^{-1}) @ C = [x, y, z]

Matrix product is associative (note the ordering change):

([sx, sy, sz] @ P) @ (P^{-1} @ C) = [x, y, z]
[sz, sy, sx] @ (P^{-1} @ C) = [x, y, z]

P is actually unitary (pure rotation without scaling, here a simple "swap") so the inverse is also the transpose:

[sz, sy, sx] @ (P^T @ C) = [x, y, z]

And we need to apply P on the right-hand side to keep the same [w, v, u] ordering between cartesian and scaled:

[sz, sy, sx] @ (P^T @ C) @ P = [x, y, z] @ P
[sz, sy, sx] @ (P^T @ C @ P) = [z, y, x]
[sz, sy, sx] @ C' = [z, y, x]

with C' = P^T @ C @ P

You could easily verify that P = P^T = [[0, 0, 1], [0, 1, 0] , [1, 0, 0]] and therefore that C' is what I claimed it is (and upper-triangular).


Cheers,
Gaël


________________________________________
From: Shi Li <sli259 at tamu.edu>
Sent: 21 December 2020 15:09
To: Gaël Donval
Subject: Re: [ase-users] MD simulation with NPT




Thank you Gaël,

I swapped the atoms positions. But for the cell transfer, swapping the first and last row won’t change it to an upper-triangular matrix. Should I also swap the first and last column?

Best,
Shi



> On Dec 19, 2020, at 8:00 PM, Gaël Donval <ggd21 at bath.ac.uk> wrote:
>
> Dear Shi,
>
> As stated in the error, your cell matrix is not upper-triangular (it is lower-triangular). An easy way to solve the problem would be swapping x and z (in atoms.positions) and swap your cell's first and last row.
>
> Best,
> Gaël
>
>
> Dear ASE users,
>
> I am trying to apply a NPT simulation on a molecular crystal cell. I used the npt.NPT function and here is part of the script:
>
> atoms = read('mol.cif')
> atoms.center()
> atoms = atoms.repeat((N,N,N))
> atoms.set_pbc(True)
> …
>
> ###NPT equilibration ####
> md = NPT(atoms, timestep= 2 * units.fs, temperature=300 * units.kB,
>         ttime=25, pfactor=75, externalstress=6.3242e-7)
>
> …
>
> This gives me an error saying that:
>
>    raise NotImplementedError("Can (so far) only operate on lists of atoms where the computational box is an upper triangular matrix.")
> NotImplementedError: Can (so far) only operate on lists of atoms where the computational box is an upper triangular matrix.
>
> Is this related with the cell I am applying? The cif file is from previous generations so the cell is not tetragonal:
>
> Cell([[14.1543, 0.0, 0.0], [-7.073274963731859, 12.255717092338744, 0.0], [0.01606649983624844, -0.009667840015844734, 17.60119001205465]])
>
> I am not sure if I did this correctly, is there a way to fix this issue? Also I am a little confuse about the externalstress for the NPT function, what should I use for 1bar?
>
> Thank you for your help.
>
> Shi