Hello, I am trying to do the following via scripting:

1. setup an Abaqus model, extracting the matrices
2. computing in Python the eigenmodes

     3. writing them to an odb

I have the following scripts, and I would like to know where the problem is, as I am getting different eigenmodes with respect to Abaqus jobs. Any input would be appreciated.

from abaqus import mdb, session
from abaqusConstants import *
import regionToolset
import mesh
import time


import numpy as np


model_name = 'BlockModel'
part_name = 'Block'
material_name = 'Steel'
section_name = 'BlockSection'
assembly_name = 'Assembly'
instance_name = part_name + '-1'


L = 1.0
H = 0.2
W = 0.2


seed_size = 0.1


E = 210e9
nu = 0.3
rho = 7850


cae_file = 'block_model.cae'
input_file = 'block_model.inp'



if model_name in mdb.models:
    del mdb.models[model_name]


model = mdb.Model(name=model_name)
s = model.ConstrainedSketch(name='__profile__', sheetSize=2.0)
s.rectangle(point1=(0.0, 0.0), point2=(L, H))


part = model.Part(name=part_name, dimensionality=THREE_D, type=DEFORMABLE_BODY)
part.BaseSolidExtrude(sketch=s, depth=W)


material = model.Material(name=material_name)
material.Elastic(table=((E, nu),))
material.Density(table=((rho,),))


model.HomogeneousSolidSection(
    name=section_name,
    material=material_name,
    thickness=None
)


cells = part.cells[:]
region = regionToolset.Region(cells=cells)
part.SectionAssignment(region=region, sectionName=section_name)


assembly = model.rootAssembly
instance = assembly.Instance(name=instance_name, part=part, dependent=ON)


faces = instance.faces
leftFaces = faces.getByBoundingBox(xMin=-1e-6, xMax=1e-6)


if len(leftFaces) == 0:
    leftFaces = [f for f in faces if abs(f.pointOn[0]) < 1e-6]


assembly.Set(name='LeftFace', faces=leftFaces)
assembly.Set(name='WholeInstance', nodes=instance.nodes)


part.seedPart(size=seed_size, deviationFactor=0.1, minSizeFactor=0.1)


elemType1 = mesh.ElemType(elemCode=C3D8, elemLibrary=STANDARD)
part.setElementType(regions=(part.cells[:],), elemTypes=(elemType1,))
part.generateMesh()


assembly.regenerate()


leftSet = assembly.sets['LeftFace']
bc_labels = [node.label for node in leftSet.nodes]
np.savetxt('leftSet.txt', bc_labels)
# model.EncastreBC(name='BC_Clamped', createStepName='Initial', region=leftSet)  # do not apply this BC here and lock it afterwards


model.FrequencyStep(name='ExtractMatrices', previous='Initial', numEigen=1)


mdb.saveAs(pathName=cae_file)
mdb.Job(name='BlockMatrixJob', model=model_name).writeInput(consistencyChecking=OFF)


input_file = 'BlockMatrixJob.inp'
with open(input_file, 'r') as f:
    lines = f.readlines()


target_line = "*Step, name=ExtractMatrices, nlgeom=NO, perturbation"
for i, line in enumerate(lines):
    if target_line in line:
        lines[i] = "*Step, name=ExtractMatrices\\n"


target_line = '*Frequency, eigensolver=Lanczos, sim, acoustic coupling=on, normalization=mass'
for i, line in enumerate(lines):
    if target_line in line:
        lines[i] = "*MATRIX GENERATE, STIFFNESS, MASS\\n"


target_line = "1, , , , ,"
for i, line in enumerate(lines):
    if target_line in line:
        lines[i] = "*MATRIX OUTPUT, STIFFNESS, MASS, FORMAT=MATRIX INPUT\\n"


with open(input_file, 'w') as f:
    f.writelines(lines)


time.sleep(5)


mdb.JobFromInputFile(name='BlockMatrixJob', inputFileName=input_file)
mdb.jobs['BlockMatrixJob'].submit()
mdb.jobs['BlockMatrixJob'].waitForCompletion()

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import eigsh


bc_labels = np.loadtxt('leftSet.txt')
bc_labels = bc_labels.astype(int)-1
bc_dofs = []
n_dofs_per_node = 3
for label in bc_labels:
    for dof in range(n_dofs_per_node):
        bc_dofs.append(label * n_dofs_per_node + dof)


data = np.loadtxt('BlockMatrixJob_STIF1.mtx', delimiter=',')


node_i = data[:,0].astype(int) - 1
dof_i  = data[:,1].astype(int) - 1
node_j = data[:,2].astype(int) - 1
dof_j  = data[:,3].astype(int) - 1
vals   = data[:,4]


n_dofs_per_node = 3
rows = node_i * n_dofs_per_node + dof_i
cols = node_j * n_dofs_per_node + dof_j


K = coo_matrix((vals, (rows, cols))).tocsr()
K = K.todense()
K = 1/2 * (K + K.T)


data = np.loadtxt('BlockMatrixJob_MASS1.mtx', delimiter=',')


node_i = data[:,0].astype(int) - 1
dof_i  = data[:,1].astype(int) - 1
node_j = data[:,2].astype(int) - 1
dof_j  = data[:,3].astype(int) - 1
vals   = data[:,4]


n_dofs_per_node = 3
rows = node_i * n_dofs_per_node + dof_i
cols = node_j * n_dofs_per_node + dof_j


M = coo_matrix((vals, (rows, cols))).tocsr()
M = M.todense()


not_bc_dofs = np.setdiff1d(np.arange(K.shape[0]), bc_dofs)


Mr = M[np.ix_(not_bc_dofs, not_bc_dofs)]
Kr = K[np.ix_(not_bc_dofs, not_bc_dofs)]


vals, vecs = eigsh(Kr, 20, M=Mr)
freqs = np.sqrt(vals) / (2 * np.pi)


true_vecs = np.zeros((K.shape[0], vecs.shape[1]))
true_vecs[not_bc_dofs, :] = vecs


np.savetxt('eigenmodes.txt', true_vecs)
np.savetxt('eigenfreqs.txt', freqs)