ElasticTensor bug when use_pressure=False

[Leimeroth]

Summary Elastic Tensor jobs work using pressure, but not using energy my strructures. When the "use_pressure" = False option is set they throw:

_fit_coeffs_with_energies(strain, energy, volume, rotations, max_polynomial_order, fit_first_order)
    150     coeff = np.triu(np.ones((6,6))).flatten()
    151     # Multiply upper triangle with upper triangle coeffs (v.s.)
--> 152     coeff[coeff!=0] *= reg.coef_[:21]*eV_div_A3_to_GPa
    153     coeff = coeff.reshape(6,6)
    154     coeff = 0.5*(coeff+coeff.T)

ValueError: operands could not be broadcast together with shapes (21,) (18,) (21,) 

Steps to Reproduce

Run a ElasticTensor calculation with the attached structure (lammps format, original atom types are Cu and Zr, but that shouldn't matter I guess) using the use_pressure=False option. structure.txt

[samwaseda]

I can't reproduce the error (and also I don't really understand how it could possibly happen). Can you maybe give more lines?

[Leimeroth]

This is what I did to test this:

ref = pr.create_job("Lammps", "ref")
ref.structure = ref_struct
ref.potential = ref.list_potentials()[0]

job = ref.create_job("ElasticTensor", "test")
job.input["use_pressure"] = False
job.run()

job2 = ref.create_job("ElasticTensor", "test2")
job2.run()

The first one fails with the error from above, while the 2nd works.

[Leimeroth]

This problem just came across my way again, this time with simple Zr. Than I ran it again, first in another project, than in the same one (using delete_existing_job=True in pr.create.job.ElasticTensor) and both times it worked. Is it possible that there is some data race / function not falways completely executed / whatever leading to random failures?

EDIT: comparing the elastic tensors produced with both jobs also give heavily different results. F.e the c66 component differs by a factor of 2. or 60 GPa. Could you please have a look at this @samwaseda

What I did:

pr = Project("ElasticTensortTest/")
lmp = pr.create.job.Lammps("ref")
lmp.structure = pr.create.structure.ase.bulk("Zr")
lmp.potential = '2015--Wilson-S-R--Ni-Zr--LAMMPS--ipr1'

et_stress = pr.create.job.ElasticTensor("ZrStress")
et_stress.ref_job = lmp
et_stress.run()

et_energy = pr.create.job.ElasticTensor("ZrEnergy")
et_energy.input["use_pressure"] = False
et_energy.ref_job = lmp
et_energy.run()

Full Backtrace:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [16], in <module>
      3 et_energy.input["use_pressure"] = False
      4 et_energy.ref_job = lmp
----> 5 et_energy.run()

File ~/git_projects/pyiron_base/pyiron_base/generic/util.py:213, in Deprecator.__deprecate_argument.<locals>.decorated(*args, **kwargs)
    203     if kw in self.arguments:
    204         warnings.warn(
    205             message_format.format(
    206                 "{}.{}({}={})".format(
   (...)
    211             stacklevel=2,
    212         )
--> 213 return function(*args, **kwargs)

File ~/git_projects/pyiron_base/pyiron_base/job/generic.py:689, in GenericJob.run(self, delete_existing_job, repair, debug, run_mode, run_again)
    687     self._run_if_repair()
    688 elif status == "initialized":
--> 689     self._run_if_new(debug=debug)
    690 elif status == "created":
    691     self._run_if_created()

File ~/git_projects/pyiron_base/pyiron_base/job/generic.py:1399, in GenericJob._run_if_new(self, debug)
   1397 else:
   1398     self.save()
-> 1399     self.run()

File ~/git_projects/pyiron_base/pyiron_base/generic/util.py:213, in Deprecator.__deprecate_argument.<locals>.decorated(*args, **kwargs)
    203     if kw in self.arguments:
    204         warnings.warn(
    205             message_format.format(
    206                 "{}.{}({}={})".format(
   (...)
    211             stacklevel=2,
    212         )
--> 213 return function(*args, **kwargs)

File ~/git_projects/pyiron_base/pyiron_base/job/generic.py:691, in GenericJob.run(self, delete_existing_job, repair, debug, run_mode, run_again)
    689     self._run_if_new(debug=debug)
    690 elif status == "created":
--> 691     self._run_if_created()
    692 elif status == "submitted":
    693     self._run_if_submitted()

File ~/git_projects/pyiron_base/pyiron_base/job/generic.py:1418, in GenericJob._run_if_created(self)
   1416     self.run_if_manually(_manually_print=False)
   1417 elif self.server.run_mode.modal:
-> 1418     self.run_static()
   1419 elif (
   1420     self.server.run_mode.non_modal
   1421     or self.server.run_mode.thread
   1422     or self.server.run_mode.worker
   1423 ):
   1424     self.run_if_non_modal()

File ~/git_projects/pyiron_base/pyiron_base/master/parallel.py:661, in ParallelMaster.run_static(self)
    657     self._run_if_master_non_modal_child_non_modal(job)
    658 elif (self.server.run_mode.modal and job.server.run_mode.modal) or (
    659     self.server.run_mode.interactive and job.server.run_mode.interactive
    660 ):
--> 661     self._run_if_master_modal_child_modal(job)
    662 elif self.server.run_mode.modal and job.server.run_mode.non_modal:
    663     self._run_if_master_modal_child_non_modal(job)

File ~/git_projects/pyiron_base/pyiron_base/master/parallel.py:584, in ParallelMaster._run_if_master_modal_child_modal(self, job)
    582 if self.is_finished():
    583     self.status.collect = True
--> 584     self.run()
    585 elif self.status.busy:
    586     self.status.refresh = True

File ~/git_projects/pyiron_base/pyiron_base/generic/util.py:213, in Deprecator.__deprecate_argument.<locals>.decorated(*args, **kwargs)
    203     if kw in self.arguments:
    204         warnings.warn(
    205             message_format.format(
    206                 "{}.{}({}={})".format(
   (...)
    211             stacklevel=2,
    212         )
--> 213 return function(*args, **kwargs)

File ~/git_projects/pyiron_base/pyiron_base/job/generic.py:697, in GenericJob.run(self, delete_existing_job, repair, debug, run_mode, run_again)
    695     self._run_if_running()
    696 elif status == "collect":
--> 697     self._run_if_collect()
    698 elif status == "suspend":
    699     self._run_if_suspended()

File ~/git_projects/pyiron_base/pyiron_base/master/parallel.py:470, in ParallelMaster._run_if_collect(self)
    462 """
    463 Internal helper function the run if collect function is called when the job status is 'collect'. It collects
    464 the simulation output using the standardized functions collect_output() and collect_logfiles(). Afterwards the
    465 status is set to 'finished'.
    466 """
    467 self._logger.info(
    468     "{}, status: {}, finished".format(self.job_info_str, self.status)
    469 )
--> 470 self.collect_output()
    472 job_id = self.get_job_id()
    473 db_dict = {}

File ~/git_projects/pyiron_atomistics/pyiron_atomistics/atomistics/master/elastic.py:561, in ElasticTensor.collect_output(self)
    559 if not self.input["use_pressure"]:
    560     output_data["pressures"] = []
--> 561 elastic_tensor, score = calc_elastic_tensor(
    562     strain=self.input["strain_matrices"],
    563     stress=-np.array(output_data["pressures"]),
    564     rotations=self.input["rotations"],
    565     energy=energy,
    566     volume=output_data["volume"],
    567     return_score=True,
    568     fit_first_order=self.input["fit_first_order"],
    569     max_polynomial_order=self.input["polynomial_order"],
    570 )
    571 output_data["fit_score"] = score
    572 output_data["elastic_tensor"] = elastic_tensor

File ~/git_projects/pyiron_atomistics/pyiron_atomistics/atomistics/master/elastic.py:299, in calc_elastic_tensor(strain, stress, energy, volume, rotations, return_score, max_polynomial_order, fit_first_order)
    286     coeff, score = _fit_coeffs_with_stress(
    287         strain=strain,
    288         stress=stress,
   (...)
    291         fit_first_order=fit_first_order,
    292     )
    293 elif (
    294     energy is not None
    295     and volume is not None
    296     and len(energy) == len(strain)
    297     and len(energy) == len(volume)
    298 ):
--> 299     coeff, score = _fit_coeffs_with_energies(
    300         strain=strain,
    301         energy=energy,
    302         volume=volume,
    303         rotations=rotations,
    304         max_polynomial_order=max_polynomial_order,
    305         fit_first_order=fit_first_order,
    306     )
    307 else:
    308     raise ValueError(
    309         "Provide either strain and stress, or strain, energy "
    310         + "and volume of same length."
    311     )

File ~/git_projects/pyiron_atomistics/pyiron_atomistics/atomistics/master/elastic.py:165, in _fit_coeffs_with_energies(strain, energy, volume, rotations, max_polynomial_order, fit_first_order)
    163 coeff = np.triu(np.ones((6, 6))).flatten()
    164 # Multiply upper triangle with upper triangle coeffs (v.s.)
--> 165 coeff[coeff != 0] *= reg.coef_[:21] * eV_div_A3_to_GPa
    166 coeff = coeff.reshape(6, 6)
    167 coeff = 0.5 * (coeff + coeff.T)

ValueError: operands could not be broadcast together with shapes (21,) (20,) (21,) 

[samwaseda]

AHA! I guess now I know where it comes from. Thanks for letting me know!

[samwaseda]

I kind of updated it. Can you give it a try?

[Leimeroth]

Yes give me a second. But since the errors seems to occur randomly I am not sure how to test it reliably, Could you explain a bit what this changes and if this also adresses my edit?

[Leimeroth]

So I tried a few times and it didn't throw the error. But the difference between elastic tensor calculated using energy and the one calculated using pressure still seems too large

[Leimeroth]

Also (c11-c12) / 2 != c66 with a non negligible difference I would say, which should be the case for hexagonal systems, if that is of any help

[samwaseda]

ooook let me see

[samwaseda]

So I think this time it should work. Can you give it a try?

[samwaseda]

Here's what was the initial problem: Imagine an input data of n-points. In order to fit the elastic constants, I first created the strain dataset by:

[[[s_1_1 * s_1_1, ..., s_1_1 * s_1_6],
...
[s_1_6 * s_1_1, ..., s_1_6 * s_1_6]],
.........
[[s_n_1 * s_n_1, ..., s_n_1 * s_n_6],
...
[s_n_6 * s_n_1, ..., s_n_6 * s_n_6]]]

where s_i_j is the j-th component of the strain vector (Voigt notation) of the i-th data point. Since the tensor is symmetric, it would only confuse the linear regression if the whole tensor is flattened and inserted as the input. In order to avoid this, I multiplied each entry by an upper triangular matrix mask, i.e. with a matrix which contains 1 in the upper triangle and 0 in the lower triangle, so that the dataset became like this:

[[[s_1_1 * s_1_1, ..., s_1_1 * s_1_6],
...
[0, ..., s_1_6 * s_1_6]],
.........
[[s_n_1 * s_n_1, ..., s_n_1 * s_n_6],
...
[0, ..., s_n_6 * s_n_6]]]

Then I flattened the array in the last two axes and removed the columns, where the sum was equal to zero, following the logic that the sum of 0 is always 0. Now the problem comes from the fact that due to symmetries, it was possible to have columns, where the sum turned to be 0 and they were removed as well, and therefore came the shape problem as you witnessed in the error messages.

[Leimeroth]

The problem that elastic tensor calculated using stress differs massively from the one calculated using the energy is still there. To make sure I also tried with another potential.

This is what I did. Differences are > 100GPa for c11

pr = Project("Testing/ElasticFix")
pr.remove_jobs(silently=True, recursive=True)

pot = '2016--Borovikov-V--Cu-Zr--LAMMPS--ipr1'

relax = pr.create.job.Lammps("Relax")
relax.structure = pr.create.structure.ase.bulk("Zr")
relax.potential = pot
relax.calc_minimize(pressure=np.zeros(6))
relax.run()

lmp = pr.create.job.Lammps("ref")
lmp.structure = relax.get_structure()
lmp.potential = pot

et_stress = pr.create.job.ElasticTensor("ZrStress")
et_stress.ref_job = lmp
et_stress.run()

et_energy = pr.create.job.ElasticTensor("ZrEnergy")
et_energy.input["use_pressure"] = False
et_energy.ref_job = lmp
et_energy.run()

ete = et_energy["output/elastic_tensor"]
ets = et_stress["output/elastic_tensor"]
ets - ete

[samwaseda]

Hm that's weird, because in my case the difference is less than 10 GPa... Let me see if I can find it out.

[samwaseda]

It looks like if orthorhombic=True is set in the structure, it gives the correct values. I just don't know why.

[Leimeroth]

Could this be related to the problems with different voigt stress definitions (i.e. if correct order is xy, xz, yz or yz, xz, xy) that cause trouble f.e. when parsing vasp output? Didn't have a look at this yet.

[samwaseda]

Could this be related to the problems with different voigt stress definitions (i.e. if correct order is xy, xz, yz or yz, xz, xy) that cause trouble f.e. when parsing vasp output? Didn't have a look at this yet.

That's rather unlikely because if that was the case orthorhombic shouldn't play a role.

[pmrv]

While you fix this: What's a known good way using this in the meantime? I've tried combinations now of non/orthorhombic, use_pressure=True/False with lammps and get different results and even non-symmetric tensors sometimes (with a minimized structure).

[samwaseda]

Did you do it with the changes from this PR? (Actually this discussion should continue in the PR)