<jupyter_start><jupyter_text>Ni - X computations part 1 Setting up the Onsager Calculator<jupyter_code>import sys
sys.path.append("../")
import numpy as np
from Onsager_calc_db import *
import onsager.crystal as crystal
from states import *
from stars import *
from vector_stars import *
import time
import pickle

# Note - a deprecation warning may be printed, but it related to a cluster expansion code built into the
# onsager code (Trinkle-2017, Trinkle-2018) and does not affect this extension of onsager for dumbbell
# mediated diffusion.
# make FCC lattice - all length units we use are in nanometers.
latt = np.array([[0., 0.5, 0.5], [0.5, 0., 0.5], [0.5, 0.5, 0.]]) * 0.352
Ni = crystal.Crystal(latt, [[np.array([0., 0., 0.])]], ["Ni"])

# Next, we set up the dumbbell containers.
# The "pdbcontainer" and "mdbcontainer" objects (see states.py) contain all information regarding pure and
# mixed dumbbells, including: The possible orientations (pdbcontainer.iorlist), the symmetry grouping
# of the dumbbells (pdbcontainer.symorlist), the group operations between them, and functions to generate
# pure and mixed dumbbell jump networks.

# We set up orientation vectors along <100> directions.
# Note that the orientation vector is a nominal vector. It helps in symmetry analysis and
# identificiation of jump types, but atomic displacements are only considered site-to-site in the
# transport coefficients. This simplification is allowed by the invariance principle of mass
# transport (Trinkle, 2018).
# To keep things somewhat physically relevant, we choose this orientation vector length to be the same
# as the host atomic diameter (0.326 nm for Ni).
o = np.array([1.,0.,0.])*0.163*2
famp0 = [o.copy()]  # multiple orientation families can be set up if desired, but we need only one.
family = [famp0]
pdbcontainer_Ni = dbStates(Ni, 0, family)
mdbcontainer_Ni = mStates(Ni, 0, family)

# Calculate the omega0 and omega2 jump networks.
# The inputs to both functions are as follows:
# The first input - the cutoff site-to-site distance for dumbbell jumps to be considered. We set this up to
# be the nearest neighbor distance.
# The next two inputs are cutoff distances for solute-solvent atoms and solvent-solvent atoms, such that
# if two atoms come within that distance of each other, then they are considered to collide (see collision.py)
jset0, jset2 = pdbcontainer_Ni.jumpnetwork(0.25, 0.01, 0.01), mdbcontainer_Ni.jumpnetwork(0.25, 0.01, 0.01)
print(Ni)
# Modify jnet0
jnet0 = jset0[0]
# These contain jumps in the form of jump objects

jnet0_indexed = jset0[1] 
# These contain jumps in the form of ((i, j), dx) where "i" is the index assigned to the initial dumbbell state
# and j is the index assigned to the final dumbbell state.

# The 90-degree roto-translation jump is the jump that has the shortest total atomic dispalcements.
# Let's try to sort the jumps accordingly.

# We first get rid of the rotation jumps
z = np.zeros(3)
indices = []
for jt, jlist in enumerate(jnet0):
    
    # We'll skip on-site rotations for now and add them later
    if np.allclose(jnet0_indexed[jt][0][1], z):
        continue
    indices.append(jt)
    
def sortkey(entry):
    jmp = jnet0[entry][0]
    # Get the initial orientation vector.
    or1 = pdbcontainer_Ni.iorlist[jmp.state1.iorind][1]
    
    # Get the Final orientation vector.
    or2 = pdbcontainer_Ni.iorlist[jmp.state2.iorind][1]
    
    # Get the site-to-site displacements
    dx = disp(pdbcontainer_Ni, jmp.state1, jmp.state2)
    
    # Get the nominal displacements along orientation vectors, which are only used to sort the jumps.
    # For the jump object, the quantity "c1" says whether the atom at the head (c1=1) or tail(c1=-1)
    # of the orientation vector is executing the jump, while the quantity "c2" says whether the jumping atom
    # ends up at the head (c2=1) or the tail (c2=-1) of the final orientation vector (see representations.py).
    dx1 = np.linalg.norm(jmp.c1*or1/2.)
    dx2 = np.linalg.norm(dx + jmp.c2*or2/2. - jmp.c1*or1/2.)
    dx3 = np.linalg.norm(-jmp.c2*or2/2.)
    return dx1+dx2+dx3  # return the total displacement of the atoms.

# Sort the indices according to total displacements.
ind_sort = sorted(indices, key=sortkey)
#keep only the required index
ind_sort = [ind_sort[0]]
# Add in the rotations
for jt, jlist in enumerate(jnet0):
    if np.allclose(jnet0_indexed[jt][0][1], z):
        ind_sort.append(jt)
ind_sort
# Let's check if we got the correct jump

# For a dumbbell, the (i, or) indices (iorInd) correspond to (basis site, orientation vector) pairs.
# The corresponding values can be found at onsagercalculator.pdbcontainer.iorlist, which we print out
# in the next cell.

print(jnet0[ind_sort[0]][0])
print(jnet0[ind_sort[1]][0])
# Indices are assigned to dumbbell (basis index, orientation) pairs based on their order in the following list.

pdbcontainer_Ni.iorlist
# Remake the omega0 jump networks with the jumps that we require.
jset0new = ([jnet0[i] for i in ind_sort], [jnet0_indexed[i] for i in ind_sort])
# Now, we modify the mixed dumbbell jumpnetwork to also give the lowest displacement jump
# Modify jnet2
jnet2 = jset2[0]
jnet2_indexed = jset2[1]
# Let's try to sort the jumps according to closest distance
# we don't want the rotational jumps as before.
z = np.zeros(3)
indices2 = []
for jt, jlist in enumerate(jnet2):
    if np.allclose(jnet2_indexed[jt][0][1], z):
        continue
    indices2.append(jt)    
print(indices2)

def sortkey2(entry):
    jmp = jnet2[entry][0]
    or1 = mdbcontainer_Ni.iorlist[jmp.state1.db.iorind][1]
    or2 = mdbcontainer_Ni.iorlist[jmp.state2.db.iorind][1]
    dx = disp(mdbcontainer_Ni, jmp.state1, jmp.state2)
    # c1 and c2 are always +1 for mixed dumbbell jumps.
    dx1 = np.linalg.norm(jmp.c1*or1/2.)
    dx2 = np.linalg.norm(dx + jmp.c2*or2/2. - jmp.c1*or1/2.)
    dx3 = np.linalg.norm(-jmp.c2*or2/2.)
    return dx1+dx2+dx3

ind_sort2 = sorted(indices2, key=sortkey2)
ind_sort2 = [ind_sort2[0]]
# Add in the rotations
for jt, jlist in enumerate(jnet2):
    if np.allclose(jnet2_indexed[jt][0][1], z):
        ind_sort2.append(jt)
print(ind_sort2)
# check if we have the correct type of jump
print(jnet2[ind_sort2[0]][0])
print(jnet2[ind_sort2[1]][0])
for tup in mdbcontainer_Ni.iorlist:
    print(tup)
# Collect the jumps that we want
jset2new = ([jnet2[i] for i in ind_sort2], [jnet2_indexed[i] for i in ind_sort2])
# Next, make an initial onsager calculator
# All possible omega4-3 jumps will first be found, the cutoff distances of which are indicated as follows:
# 0.25 : cutoff for site-to-site atomic displacement during omega4-3 jumps.
# 0.01 : distance of closest approach during omega4-3 jumps. Input twice for solute-solute and solute-solvent
# jumps.
# Nthermo - the thermodynamic shell - one nearest neighbor in the present case.
start = time.time()
onsagercalculator = dumbbellMediated(pdbcontainer_Ni, mdbcontainer_Ni, jset0new, jset2new, 0.25,
                                     0.01, 0.01, 0.01, NGFmax=4, Nthermo=1)
print("onsager calculator initiation time = {}".format(time.time() - start))
# Next, we get the omega43 jumps computed internally in the onsager calculator,
# and extract the shortest displacement jumps in the same manner as the omega0 and omega2 jumps.
jnet43 = onsagercalculator.jnet43
jnet43_indexed = onsagercalculator.jnet43_indexed
# Let's try to sort the jumps according to closest distance
z = np.zeros(3)
indices43 = []
for jt, jlist in enumerate(jnet43):
    if np.allclose(jnet43_indexed[jt][0][1], z):
        continue
    indices43.append(jt)    
print(indices43)

def sortkey43(entry):
    jmp = jnet43[entry][0] # This is an omega4 jump
    if not jmp.c2 == -1:
        print(c2)
    or1 = pdbcontainer_Ni.iorlist[jmp.state1.db.iorind][1]
    or2 = mdbcontainer_Ni.iorlist[jmp.state2.db.iorind][1]
    dx = disp4(pdbcontainer_Ni, mdbcontainer_Ni, jmp.state1, jmp.state2)
    # remember that c2 is -1 for an omega4 jump
    dx1 = np.linalg.norm(jmp.c1*or1/2.)
    dx2 = np.linalg.norm(dx - or2/2. - jmp.c1*or1/2.)
    dx3 = np.linalg.norm(jmp.c2*or2/2.)
    return dx1+dx2+dx3

ind_sort43 = sorted(indices43, key=sortkey43)[:1]
print(ind_sort43)
# check if we have the correct jump

# Next, we take a look at our 90-degree roto-translational omega4-omega3 jumps
# Since they have the same transition state energy, they are all grouped together in
# the same list.

# The Initial state is the complex state, and the Final state is the mixed dumbbell state involved
# in the omega4 jump (the reverse is the omega3 jump.

# the list into which these states are indexed are printed in the following cells


print(jnet43[ind_sort43[0]][0])
print(jnet43_indexed[ind_sort43[0]][0])
# The pure dumbbell state corresponding to the (i, or) index of the Initial state is in the same index
# in the following list
pdbcontainer_Ni.iorlist
# The mixed dumbbell state corresponding to the (i, or) index of the Final state is in the same index
# place in the following list
mdbcontainer_Ni.iorlist
# reconstruct the omega43 jump network to include only the 90-degree roto-translational jumps we have just
# filtered out.
onsagercalculator.regenerate43(ind_sort43)
# We then store the onsager calculator into a pickle file so that it need not be regenerated all the time.
with open('NiFe_NiCr_Onsg.pkl','wb') as fl:
    pickle.dump(onsagercalculator,fl)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09574920297555792, 'fraction_numerical': 0.03623804463336876, 'mean_word_length': 4.092532467532467, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '268061', 'license': 'MIT', 'n_tokens_mistral': 3387, 'n_tokens_neox': 3118, 'n_words': 1136, 'path': 'FCC_Calculations/.ipynb_checkpoints/Ni_X_Par1_setup-checkpoint.ipynb', 'repo_name': 'sohamch/dumbbell_functions_initial', 'size': '9410'}