fixed tbias

# LAMMPS GJF-2GJ THERMOSTAT EXAMPLE W/ PYTHON

## GJF-2GJ THERMOSTAT

This directory contains a python script to run NVT simulations using the GJF-2GJ thermostat.

The script will vary the timestep and write thermodynamic output to screen.

This script has True/False options to change how you would like to dump/write your output.

Example:

```

NP=4 #number of processors

mpirun -np $NP python gjf.py

```

## Required LAMMPS packages: MOLECULE package

## LAMMPS COMPILE MODE: SHLIB

## LAMMPS OPTIONAL INSTALL: make install-python

## Required Python packages: mpi4py 

#############################

#Atoms types - mass - charge#

#############################

#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#

variable Ar equal 1

#############

#Atom Masses#

#############

mass ${Ar}      39.903

###########################

#Pair Potentials - Tersoff#

###########################

pair_style      lj/cubic

pair_coeff      * * 0.0102701 3.42

"""Made by Charlie Sievers Ph.D. Candidate, UC Davis, Donadio Lab 2019"""

from mpi4py import MPI

from lammps import lammps

import lammps_tools as lt

import numpy as np

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

""" LAMMPS  VARIABLES """

# new file or restart

run_no = 0

# data files

infile = "argon.lmp"

restart_file = "final_restart.{}".format(run_no)

ff_file = "ff-argon.lmp"

outfile = "output.dat"

# write final_restart

write_final_restart = False

# random numbers

seed0 = 2357

seed1 = 26588

seed2 = 10669

# MD Parameters

# number of steps

nsteps = 50000

# timestep

# dt = 0.001

# starting simulation temp

temp_start = 10

# final simulation temp

temp_final = 10

# relaxation time

trel = 1

# trajectory frequency

ntraj = 0

# Ensemble 0 = GJF u, 1 = GJF v, 2 = Nose-Hoover, 3 = Langevin, 4 = BDP (Currently all NVT)

ensemble = 0

# Output Parameters

nthermo = 200

nout = int(nsteps / nthermo)  # Important

# output to screen and log file?

lammps_output = False

# Lammps Thermo

thermo = False

python_output = True

# Write output to file?

write_output = False

if write_output is True:

    data = open("{}".format(outfile), "w")

if python_output is True:

    if rank == 0:

        print("dt, temp, ke, fke, pe, fpe")

for j in range(20):

    # timestep

    dt = 0.005*(j+1)

    if lammps_output is True:

        lmp = lammps()

    else:

        lmp = lammps(cmdargs=["-screen", "none", "-log", "none"])

    lmp.command("atom_style full")

    lmp.command("units metal")

    lmp.command("processors * * *")

    lmp.command("neighbor 1 bin")

    lmp.command("boundary p p p")

    if run_no is 0:

        lmp.command("read_data {}".format(infile))

    else:

        lmp.command("read_restart final_restart".format(run_no-1))

    if thermo is True:

        lmp.command("thermo_style custom time temp  pe ke press vol  cpu")

        lmp.command("thermo {}".format(nthermo))

        lmp.command("thermo_modify flush yes")

    lmp.file("{}".format(ff_file))

    lmp.command("timestep {}".format(dt))

    # get_per_atom_compute example with dim of two and within a group

    # lmp.command("region rand block 5 20 5 20 5 20")

    # lmp.command("group rand region rand")

    # lmp.command("compute x rand property/atom x y")

    # test = get_per_atom_compute(comm, lmp, "x", 2, group="rand")

    lmp.command("compute ke all ke/atom")

    lmp.command("compute pe all pe")

    if ntraj != 0:

        lmp.command("dump 1 all dcd {} trajectory.dcd".format(ntraj))

        lmp.command("dump_modify 1 unwrap yes")

    if run_no == 0:

        lmp.command("velocity all create {} {} mom yes dist gaussian".format(temp_start, seed0))

    lmp.command("fix nve all nve")

    if ensemble == 0:

        # gjf u

        lmp.command("fix lang all langevin {} {} {} {} gjf yes halfstep yes".format(

            temp_start, temp_final, trel, seed1))

    elif ensemble == 1:

        # gjf v

        lmp.command("fix lang all langevin {} {} {} {} gjf yes".format(

            temp_start, temp_final, trel, seed1))

    elif ensemble == 2:

        # NH

        lmp.command("fix nvt all nvt temp {} {} {}".format(

            temp_start, temp_final, trel))

    elif ensemble == 3:

        # lang

        lmp.command("fix lang all langevin {} {} {} {} tally yes zero yes".format(

            temp_start, temp_final, trel, seed1))

    elif ensemble == 4:

        # BDP

        lmp.command("fix stoch all temp/csvr {} {} {} {}".format(

            temp_start, temp_final, trel, seed1))

    natoms = lmp.extract_global("natoms", 0)

    nlocal = lmp.extract_global("nlocal", 0)

    ke_sum = lt.get_per_atom_compute(comm, lmp, "ke")

    ke_2 = ke_sum**2

    pe_sum = 0

    pe_2 = 0

    temp_sum = 0

    for i in range(nout):

        nlocal = lmp.extract_global("nlocal", 0)

        lmp.command("run {} pre no post no".format(nthermo))

        temp = lmp.extract_compute("thermo_temp", 0, 0)

        ke = lt.get_per_atom_compute(comm, lmp, "ke")

        pe = lmp.extract_compute("pe", 0, 0)

        ke_sum += ke

        ke_2 += ke**2

        pe_sum += pe

        pe_2 += pe**2

        temp_sum += temp

        if python_output is True:

            if rank == 0:

                print("Time: {:.6f}, Temp: {:.6f}, KE: {:.6f}, PE: {:.6f}".format(

                    i*nthermo*dt, temp, ke.sum(), pe))

    if write_final_restart is True:

        lmp.command("write_restart {}".format(restart_file))

    if rank == 0:

        ke = ke_sum.sum() / (nout + 1)

        fke = (np.sqrt((ke_2 - ke_sum ** 2 / (nout + 1)) / (nout + 1))).sum()

        pe = pe_sum / nout

        fpe = np.sqrt((pe_2 - pe_sum ** 2 / nout) / nout)

        temp = temp_sum / nout

        if python_output is True:

            print(dt, temp, ke, fke, pe, fpe)

        if write_output is True:

            data.write("{:.6f} {:.6f} {:.6f} {:.6f} {:.6f} {:.6f}\n".format(

                dt, temp, ke, fke, pe, fpe))

            data.flush()

if write_output is True:

    data.close()

"""Made by Charlie Sievers Ph.D. Candidate, UC Davis, Donadio Lab 2019"""

from mpi4py import MPI

import numpy as np

import ctypes as ctypes

""" USEFULL LAMMPS FUNCTION """

def get_nlocal(lmp):

    nlocal = lmp.extract_global("nlocal", 0)

    return nlocal

def get_aid(lmp, group=None):

    if group is None:

        c_aid = lmp.extract_atom("id", 0)

        ptr = ctypes.cast(c_aid, ctypes.POINTER(ctypes.c_int32 * get_nlocal(lmp)))

        aid = np.frombuffer(ptr.contents, dtype=np.int32)

    else:

        try:

            c_aid = lmp.extract_variable("aid", group, 1)

            ptr = ctypes.cast(c_aid, ctypes.POINTER(ctypes.c_double * get_nlocal(lmp)))

            aid = np.frombuffer(ptr.contents, dtype=np.double)

        except ValueError:

            lmp.command("variable aid atom id")

            aid = get_aid(lmp, group)

    return aid

def get_per_atom_compute(comm, lmp, name, dim=1, dtype="double", group=None):

    laid = get_aid(lmp, group)

    nlocal = get_nlocal(lmp)

    ngroup = comm.allgather(laid)

    type = dim

    if dim > 1:

        type = 2

    for array in ngroup:

        try:

            aid = np.concatenate((aid, array))

        except UnboundLocalError:

            aid = array

    if dtype == "double":

        mem_type = ctypes.c_double

    elif dtype == "integer":

        mem_type = ctypes.c_int

    elif dtype == "bigint":

        mem_type = ctypes.c_int32

    else:

        print("{} not implemented".format(dtype))

        return

    tmp = lmp.extract_compute(name, 1, type)

    if type == 1:

        ptr = ctypes.cast(tmp, ctypes.POINTER(mem_type * nlocal))

    else:

        ptr = ctypes.cast(tmp[0], ctypes.POINTER(mem_type * nlocal * dim))

    lcompute = comm.allgather(np.frombuffer(ptr.contents).reshape((-1, dim)))

    for array in lcompute:

        try:

            compute = np.concatenate((compute, array))

        except UnboundLocalError:

            compute = array

    aid = np.expand_dims(aid, axis=1)

    compute = np.concatenate((aid, compute), axis=-1)

    compute = compute[compute[..., 0] != 0]

    compute = compute[compute[..., 0].argsort()][..., 1:]

    if dim == 1:

        compute = np.squeeze(compute, axis=-1)

    return compute

 No newline at end of file