2
$\begingroup$

How can one verify that the values (0, 1, 2, and 3) in the SMASS parameter correspond appropriately to the NVT ensemble(Nose-Hoover thermostat) in the context of an ab initio molecular dynamics simulation on a Pd(111) surface with dimensions 2x2x6 at 300K?

$\endgroup$

1 Answer 1

5
$\begingroup$

If you are on VASP, there is a script from Qijing Zheng to whom goes the credit, with this script you can estimate the value of SMASS. Just in case you can't access that, I am leaving you the script in here:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import argparse
import ase
import sys
import numpy as np
from ase.io import read, write


def nose_mass(temperature, ndof, t0, L):
    '''
    Suggested Q:

        Shuichi Nosé, J. Chem. Phys., 81, 511(1984). 

    input:

    temperaute: in unit of Kelvin
    ndof: No. of degrees of freedom
    t0: The oscillation time in fs
    L: the length of the first basis vector
    '''

    # Q in Energy * Times**2
    qtmp = (t0 * 1E-15 / np.pi / 2)**2 * \
        2 * ndof * ase.units.kB * temperature \
        * ase.units._e

    # Q in AMU * Angstrom**2
    Q = qtmp / ase.units._amu / (L * 1E-10)**2

    return Q


def cnt_dof(atoms):
    '''
    Count No. of Degrees of Freedom
    '''
    if atoms.constraints:
        from ase.constraints import FixAtoms, FixScaled, FixedPlane, FixedLine
        sflags = np.zeros((len(atoms), 3), dtype=bool)
        for constr in atoms.constraints:
            if isinstance(constr, FixScaled):
                sflags[constr.a] = constr.mask
            elif isinstance(constr, FixAtoms):
                sflags[constr.index] = [True, True, True]
            elif isinstance(constr, FixedPlane):
                mask = np.all(np.abs(np.cross(constr.dir, atoms.cell)) < 1e-5,
                              axis=1)
                if sum(mask) != 1:
                    raise RuntimeError(
                        'VASP requires that the direction of FixedPlane '
                        'constraints is parallel with one of the cell axis')
                sflags[constr.a] = mask
            elif isinstance(constr, FixedLine):
                mask = np.all(np.abs(np.cross(constr.dir, atoms.cell)) < 1e-5,
                              axis=1)
                if sum(mask) != 1:
                    raise RuntimeError(
                        'VASP requires that the direction of FixedLine '
                        'constraints is parallel with one of the cell axis')
                sflags[constr.a] = ~mask

        return np.sum(~sflags)
    else:
        return len(atoms) * 3 - 3


def parse_cml_args(cml):
    '''
    CML parser.
    '''
    arg = argparse.ArgumentParser(add_help=True)

    arg.add_argument('-i', dest='poscar', action='store', type=str,
                     default='POSCAR',
                     help='Real POSCAR to calculate the No. of Degrees of freedom')
    arg.add_argument('-u', dest='unit', action='store', type=str,
                     default='fs',
                     choices=['cm-1', 'fs'],
                     help='Default unit of input frequency')
    arg.add_argument('-T', '--temperature', dest='temperature',
                     action='store', type=float,
                     default=300,
                     help='The temperature.')
    arg.add_argument('-f', '--frequency', dest='frequency',
                     action='store', type=float,
                     default=40,
                     help='The frequency of the temperature oscillation')

    return arg.parse_args(cml)


if __name__ == '__main__':
    arg = parse_cml_args(sys.argv[1:])

    if arg.unit == 'cm-1':
        THzToCm = 33.3564095198152
        t0 = 1000 * THzToCm / arg.frequency
    else:
        t0 = arg.frequency

    pos  = read(arg.poscar)
    L    = np.linalg.norm(pos.cell, axis=1)[0]
    ndof = cnt_dof(pos)
    Q    = nose_mass(arg.temperature, ndof, t0, L)

    print("SMASS = {}".format(Q))
$\endgroup$
7
  • $\begingroup$ Thanks a lot @Zafar Mehrez but I have a doubt in this code. how to calculate the frequency of the temperature oscillation and here 40 is given as default value . $\endgroup$ Jan 23 at 9:17
  • $\begingroup$ @RubiAgrawal, as it says, it is the oscillation of temperature, you can change with -f tag in accordance to your system. If SMASS is too small then this frequency will get higher, ideally you want to be equal to the characteristic vibrational frequencies (phonon frequencies) of the real system. If you find that not suitable to use, I remember a hack to get SMASS, you can run your system using SMASS=0, IBRION=0 and TEBEG=300, and TEEND=300, NSW=1000, NBLOCK=4 after this short calculation you can find the right SMASS in the OUTCAR file $\endgroup$ Jan 23 at 14:15
  • $\begingroup$ Thank you @Jaafar Mehrez, i have one last question, Should I use different SMASS values when simulating the same system at different temperatures (cause SMASS chnage with temperature), or is it preferable to maintain the same SMASS for comparing the effects under varying temperature conditions? $\endgroup$ Jan 23 at 15:25
  • $\begingroup$ @RubiAgrawal, to give you more sense of what SMASS is, check doi.org/10.1063/1.447334, Nosé defined Q=2gK_bT/(w_0)^2. In VASP, this Q is controlled with SMASS, thus Q or SMASS controls the flucuation in temperature, I am not sure what information you could get by keeping SMASS constant for different temperatures $\endgroup$ Jan 24 at 3:28
  • $\begingroup$ @ Jaafar Mehrez I want to calculate molecular dissociation on metal surface at different temperature so for this at different temperature, So for this at different temperatures should I use the same SMASS or different ones? $\endgroup$ Jan 24 at 5:42

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .