Source code for abics.scripts.abicsRXsepT

# ab-Initio Configuration Sampling tool kit (abICS)
# Copyright (C) 2019- The University of Tokyo
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see http://www.gnu.org/licenses/.

from mpi4py import MPI
import shutil
import os, sys
import datetime
import argparse
from abics.mc_mpi import RXParams, RX_MPI_init
import numpy as np
import scipy.constants as constants


[docs]def main(): # Input parser parser = argparse.ArgumentParser( description="Reorganize abICS RXMC results by temperature" ) parser.add_argument("inputfi", help="toml input file used for abICS run") parser.add_argument( "skipsteps", nargs="?", type=int, default=0, help="number of thermalization steps to skip in energy averaging." + " Default: 0", ) args = parser.parse_args() inputfi = args.inputfi nskip = args.skipsteps rxparams = RXParams.from_toml(inputfi) nreplicas = rxparams.nreplicas comm = RX_MPI_init(rxparams) myreplica = comm.Get_rank() if myreplica == 0: if os.path.exists("Tseparate"): shutil.move("Tseparate", "Tseparate.bak.{}".format(datetime.datetime.now())) os.mkdir("Tseparate") comm.Barrier() # Separate structure files os.mkdir(os.path.join("Tseparate", str(myreplica))) Trank_hist = np.load(os.path.join(str(myreplica), "Trank_hist.npy")) os.chdir(str(myreplica)) for j in range(len(Trank_hist)): shutil.copy( "structure.{}.vasp".format(j), os.path.join(os.pardir, "Tseparate", str(Trank_hist[j])), ) # Separate energies myreplica_energies = np.load("obs_save.npy")[:, 0] for Tid in range(nreplicas): T_energies = np.where(Trank_hist == Tid, myreplica_energies, 0) T_energies_rcvbuf = np.zeros(T_energies.shape[0], "d") comm.Reduce( [T_energies, MPI.DOUBLE], [T_energies_rcvbuf, MPI.DOUBLE], op=MPI.SUM, root=Tid, ) if myreplica == Tid: np.savetxt( os.path.join(os.pardir, "Tseparate", str(Tid), "energies.dat"), T_energies_rcvbuf, ) comm.Barrier() if myreplica == 0: os.chdir(os.path.join(os.pardir, "Tseparate")) with open("energies_T.dat", "w") as fi: kTs = np.load(os.path.join(os.pardir, "kTs.npy")) Ts = kTs / constants.value(u"Boltzmann constant in eV/K") for Tid in range(nreplicas): energy_mean = np.mean( np.loadtxt(os.path.join(str(Tid), "energies.dat"))[nskip:] ) fi.write("{}\t{}\n".format(Ts[Tid], energy_mean))
if __name__ == "__main__": main()