[ase-users] writing trajectories of a subsystem

[Ask Hjorth Larsen]

Hi,

2018-03-27 22:44 GMT+02:00 Efrem Braun via ase-users
<ase-users at listserv.fysik.dtu.dk>:
> Hi,
>
> I'm using ASE to calculate diffusion coefficients of molecules within a
> framework, and for this I only need the trajectory file to report the
> positions of the adsorbate molecules themselves. The trajectory files are
> pretty large, and so cutting out the unnecessary information of the
> framework atoms would be helpful. Is there a way to do this?
>
> I tried the following, but it didn't work. Using the example of running an
> MD simulation of a box of water molecules
> (https://wiki.fysik.dtu.dk/ase/tutorials/tipnp_equil/tipnp_equil.html), I
> made the highlighted change:
>
> ####################################################################
> from ase import Atoms
> from ase.constraints import FixBondLengths
> from ase.calculators.tip3p import TIP3P, rOH, angleHOH
> from ase.md import Langevin
> import ase.units as units
> from ase.io.trajectory import Trajectory
> import numpy as np
>
> # Set up water box at 20 deg C density
> x = angleHOH * np.pi / 180 / 2
> pos = [[0, 0, 0],
>        [0, rOH * np.cos(x), rOH * np.sin(x)],
>        [0, rOH * np.cos(x), -rOH * np.sin(x)]]
> atoms = Atoms('OH2', positions=pos)
>
> vol = ((18.01528 / 6.022140857e23) / (0.9982 / 1e24))**(1 / 3.)
> atoms.set_cell((vol, vol, vol))
> atoms.center()
>
> atoms = atoms.repeat((3, 3, 3))
> atoms.set_pbc(True)
>
> # RATTLE-type constraints on O-H1, O-H2, H1-H2.
> atoms.constraints = FixBondLengths([(3 * i + j, 3 * i + (j + 1) % 3)
>                                    for i in range(3**3)
>                                    for j in [0, 1, 2]])
>
> tag = 'tip3p_27mol_equil'
> atoms.calc = TIP3P(rc=4.5)
> md = Langevin(atoms, 1 * units.fs, temperature=300 * units.kB,
>               friction=0.01, logfile=tag + '.log')
>
> traj = Trajectory(tag + '.traj', 'w', atoms[0:10])
> # I also tried the following, with the same result:
> # atoms_subsystem = atoms[0:10]
> # traj = Trajectory(tag + '.traj', 'w', atoms_subsystem)
> md.attach(traj.write, interval=1)
> md.run(4000)
> ####################################################################
>
> The trajectory that was printed did not show any of the atoms moving from
> their initial system. It seems like the "atoms" object was copied instead of
> having atoms_subsystem pointing to a reference of the original object, even
> though https://wiki.fysik.dtu.dk/ase/ase/atom.html indicates that this
> shouldn't be.
>
> Any ideas?

You need to define a custom write function.  For example something like:

atoms = ...
traj = Trajectory(...)
md = ...

def mywritefunction():
    traj.write(atoms[:10])

md.attach(mywritefunction, interval=1)

In the commented example, you make a subsystem and write that same
subsystem many times.  But onlyy the full atoms object was actually
updated, so the initial subsystem would be copied many times into the
trajectory.

Best regards
Ask

>
> Efrem Braun
>
> _______________________________________________
> ase-users mailing list
> ase-users at listserv.fysik.dtu.dk
> https://listserv.fysik.dtu.dk/mailman/listinfo/ase-users