Fitting GAP with DFT VASP training sets

[ap439111]

Hi,

I have been looking to fit GAP from VASP MD trajectories. I have successfully installed quippy and gap and ran the tutorial for EMT. I am wondering how can I format my train.xyz to train GAP from the atom resolved forces and total energy (there are no atom split energies for DFT). If I want to use soap for the multi-component system, is specifying atomic number sufficient (e.g. Z={1,8} for H and O)? Does the order matter? Can I specify Z={8,1}? I would appreciate it if someone can link me to any tutorial for documentation to use quippy for DFT-generated data sets.

Thank you.

[gabor1]

look at the xyz generated in the tutorial (it doesn’t matter whether it is EMT or DFT), we always train from total energies and forces (and virials = -volume*stress in eV). order of specifying Z doesn’t matter. you don’t even need to specify any Z, the default is to consider all species. for other GAP parameters, please see published papers - the xml file always has the full input command line (grep for "command"). For water here is a start, which of course can be tweaked. gap_fit atoms_filename=train.xyz sparse_jitter=1e-10 gap={soap cutoff=3.0 cutoff_transition_width=1.0 n_max=8 l_max=4 delta=0.1 atom_sigma=0.3 zeta=4 n_sparse=1500 normalise=T sparse_method=cur_points add_species covariance_type=dot_product : soap cutoff=6.0 cutoff_transition_width=1.0 n_max=8 l_max=4 delta=0.05 atom_sigma=0.6 zeta=4 n_sparse=1500 normalise=T sparse_method=cur_points add_species covariance_type=dot_product} e0_method=average default_sigma={0.0005 0.02 0.0 0.0} energy_parameter_name=energy force_parameter_name=force do_copy_at_file=F sparse_separate_file=T gp_file=gp_water.xml if you put --help to any part of this instead of key=value, you get a list of allowed keys. Let me know if you encounter any difficulties.

…

-- Gábor Gábor Csányi Professor of Molecular Modelling Engineering Laboratory, University of Cambridge Pembroke College Cambridge Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

On 2 Jun 2021, at 18:47, Anup Pandey ***@***.***> wrote: Hi, I have been looking to fit GAP from VASP MD trajectories. I have successfully installed quippy and gap and ran the tutorial for EMT. I am wondering how can I format my train.xyz to train GAP from the atom resolved forces and total energy (there are no atom split energies for DFT). If I want to use soap for the multi-component system, is specifying atomic number sufficient (e.g. Z={1,8} for H and O)? Does the order matter? Can I specify Z={8,1}? I would appreciate it if someone can link me to any tutorial for documentation to use quippy for DFT-generated data sets. Thank you. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

Thank you Prof Gabor.

[ap439111]

Hi Prog Gabor,

I have successfully used the package to train the model. I have a concern while training with the virials. The training completed successfully while the virials are exactly same as the actual ones (please see the attached plots). This is highly unrealistic. Also, my assumption is that the forces are also used in the training if they are included in the xyz. Is that correct?

Thank you for your support.

Anup

[gabor1]

You are correct that predicted virials will not be exact. I suspect that your script somehow is reading the original virials. The gap code fits on forces if it finds them but in your case it didn't find any forces. Perhaps the name of the force columns in the xyz file was wrong?

…

-- Gábor

On 23 Jun 2021, at 22:08, Anup Pandey ***@***.***> wrote:  Hi Prog Gabor, I have successfully used the package to train the model. I have a concern while training with the virials. The training completed successfully while the virials are exactly same as the actual ones (please see the attached plots). This is highly unrealistic. Also, my assumption is that the forces are also used in the training if they are included in the xyz. Is that correct? Thank you for your support. Anup — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

I have forces after the coordinates (attached). The script is reading fine. I have attached 5 virials grep ed from train.xyz and quippy_train.xyz.

[gabor1]

What exact command did you use to evaluate the trained model? The forces look ok in your training set.

…

-- Gábor

On 23 Jun 2021, at 22:20, Anup Pandey ***@***.***> wrote:  I have forces after the coordinates (attached). The script is reading fine. I have attached 5 virials grep ed from train.xyz and quippy_train.xyz. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[mcaroba]

He's passing "force" instead of "forces", so they're not being read.

[ap439111]

Oh, this could be the case. Thanks !

[ap439111]

I still can't figure out the error in training the virials.

[gabor1]

As a test, use the "quip" command line code to evaluate the same xyz file

…

-- Gábor

On 23 Jun 2021, at 22:20, Anup Pandey ***@***.***> wrote:  I have forces after the coordinates (attached). The script is reading fine. I have attached 5 virials grep ed from train.xyz and quippy_train.xyz. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

quip E=T F=T atoms_filename=train.xyz param_filename=gp_rhea.xml | grep AT | sed 's/AT//' > quip_train.xyz

[gabor1]

Aha, I missed that! You are explicitly lying about the force parameter name

…

-- Gábor

On 23 Jun 2021, at 22:30, Miguel A. Caro ***@***.***> wrote:  He's passing "force" instead of "forces", so they're not being read. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[gabor1]

Right! So you are not asking to compute the virial. Hence the old virial is left in the file. Add " V " (you don't need the "=T" ) and you'll get the GAP virial

…

-- Gábor

On 23 Jun 2021, at 22:30, Anup Pandey ***@***.***> wrote:  quip E=T F=T atoms_filename=train.xyz param_filename=gp_rhea.xml | grep AT | sed 's/AT//' > quip_train.xyz — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

So, the following should be okay?

quip E F V atoms_filename=train.xyz param_filename=gp_rhea.xml | grep AT | sed 's/AT//' > quip_train.xyz

[gabor1]

Yes

…

-- Gábor

On 23 Jun 2021, at 22:39, Anup Pandey ***@***.***> wrote:  So, the following should be okay? quip E F V atoms_filename=train.xyz param_filename=gp_rhea.xml | grep AT | sed 's/AT//' > quip_train.xyz — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

Thank you !

[bernstei]

This is why I prefer to name the reference quantities something other than plain "energy", "forces", etc, (e.g. "DFT_energy"), and then get the fitting code to use those fields when fitting.

[ap439111]

That would be the right way of doing it.

I have another quick question. Is there a way to retrain the trained model?

Thanks !

[gabor1]

Not sure what you mean by retrain. we solve the linear problem with a direct method, not an iterative method, so there is no sense in which a training would continue from the result of a previous one. one could perhaps save come computation if the kernel values were available from a previous training, but we have not implemented such a workflow.

[ap439111]

Thank you.

While training from energy and forces with 940 structures I ran into memory issue. I am not sure how is it going to affect the accuracy if the training is carried out using only energy versus energy+forces. I want to eventually carry out MD using lammps, is the training with only energy sufficient ?

[gabor1]

No, you need forces. why do you think you need so many structures? how big are your structures? we typically use iterative training, start with small number, and increase if needed.

…

-- Gábor Gábor Csányi Professor of Molecular Modelling Engineering Laboratory, University of Cambridge Pembroke College Cambridge Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

On 25 Jun 2021, at 18:07, Anup Pandey ***@***.***> wrote: Thank you. While training from energy and forces with 940 structures I ran into memory issue. I am not sure how is it going to affect the accuracy if the training is carried out using only energy versus energy+forces. I want to eventually carry out MD using lammps, is the training with only energy sufficient ? — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

I want to build a model trained on configurations that are sampled through various temperatures (say 600K-1000K, at 100K steps), which would increase the config size. Each structure is 64 atoms (lattice ~ 6, 13, 13 Arng.). May be I should take few structures from each temp.

Thank you.

[gabor1]

why do you want to take so small steps in temp? I would just take 600K (or even lower) and 1000K. think about it, the kinds of configurations that are easily accessible at 600K and 700K are not THAT different!! There is typically *excellent* extrapolation in temperature. Instead of small steps in temperature, make sure you cover different volumes - the extrapolation in density is much worse. I am not a fan of taking long ab initio MDs as data. better to start with something small and using GAP md to sample where it wants to go, and fold back configs from that. (or maybe were you referring to doing this?).

…

-- Gábor Gábor Csányi Professor of Molecular Modelling Engineering Laboratory, University of Cambridge Pembroke College Cambridge Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

On 25 Jun 2021, at 18:33, Anup Pandey ***@***.***> wrote: I want to build a model trained on configurations that are sampled through various temperatures (say 600K-1000K, at 100K steps), which would increase the config size. Each structure is 64 atoms (lattice ~ 6, 13, 13 Arng.). May be I should take few structures from each temp. Thank you. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

That temperature steps are just an example. I am actually sampling at larger temp steps. yeah, I am looking to sample through GAP MD.

Thank you.

[gabor1]

OK, so if you are sampling through GAP MD, then you are doing the right thing. gallop through temperature, pressure, whatever else you can think of. but don’t add back too many things at any one step, because you will not need it!

…

-- Gábor Gábor Csányi Professor of Molecular Modelling Engineering Laboratory, University of Cambridge Pembroke College Cambridge Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

On 25 Jun 2021, at 18:47, Anup Pandey ***@***.***> wrote: That temperature steps are just an example. I am actually sampling at larger temp steps. yeah, I am looking to sample through GAP MD. Thank you. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

Thanks. I was skimming through your GST paper and realized that there were ~3000 configs in total. Are the models separately trained for melt, equilibrate, quench, Crystal etc. ?

[gabor1]

No, there was one model trained, it was one of the largest models we have trained (what matters is the total number of atoms in the training set). We have a 1.5TB memory node to train these large models on. LAter this year we will be releasing an MPI-parallelised version of the training code, which will mean that you don’t need all this memory on a single node. How much training data you need depends on the application. the GST was particularly challenging because the three components can occur virtually in any configuration within the local environment. Can I ask what your system is?

…

-- Gábor Gábor Csányi Professor of Molecular Modelling Engineering Laboratory, University of Cambridge Pembroke College Cambridge Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

On 25 Jun 2021, at 21:14, Anup Pandey ***@***.***> wrote: Thanks. I was skimming through your GST paper and realized that there were ~3000 configs in total. Are the models separately trained for melt, equilibrate, quench, Crystal etc. ? — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

Wow. I don't have access to such a large memory node. My system is high-entropy alloy and as a test case I wanna train the model on MoNbTaW. If possible, use the force field with lammps to study some mechanical properties. Any suggestion is highly appreciated. Btw for this system neural network does really bad.

[gabor1]

I predict that with the current version of SOAP, this will be difficult. Precisely because you have a nearly random alloy, so many many possibilities within the local environment, so you will indeed need a large number of training configs. The other problem is that with 4 components, the SOAP vectors (the descriptor) is large, and this exacerbates the problem. We are working hard on solving these problems, but none of these is immediate, I would say ~ 3-6 months. 1) MPI parallelisation of gap_fit, this will help with the memory 2) we have a version of SOAP that is linear scaling with the number of elements, so 4 elements will be only 4 times more expensive memory and speed wise, not 16 times, as one element. In the meantime, I have two suggestions for you. A) create a 2body+soap model, where the 2b is ~5A long, the SOAP is *nearest*neighbour* so about 3A cutoff. It will be interesting to see what accuracy you can squeeze out of such a split short range model. B) do what Byggmastar did for a 5-component system: use a 2b+3b GAP model, no SOAP. they found this to be reasonably accurate: https://arxiv.org/abs/2106.03369 In fact 4 of their 5 elements is the same as yours. just ask them for their model!

…

-- Gábor Gábor Csányi Professor of Molecular Modelling Engineering Laboratory, University of Cambridge Pembroke College Cambridge Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

On 25 Jun 2021, at 21:25, Anup Pandey ***@***.***> wrote: Wow. I don't have access to such a large memory node. My system is high-entropy alloy and as a test case I wanna train the model on MoNbTaW. If possible, use the force field with lammps to study some mechanical properties. Any suggestion is highly appreciated. Btw for this system neural network does really bad. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

Thank you.

[ap439111]

Is anyone aware of following error while using GAP with LAMMPS? I built LAMMPS with user-quip package following the instructions in a README file (/lib/quip/README) and there was no issue with the built. I tested the exec in the LAMMPS quip tutorial.

SYSTEM ABORT: quip_lammps_wrapper: quip_potential not initialised

Thank you.

[gabor1]

Is that the full output? If there is more, can you paste? best if you run serial, to make sure MPI doesn’t swallow some output lines.

…

-- Gábor Gábor Csányi Professor of Molecular Modelling Engineering Laboratory, University of Cambridge Pembroke College Cambridge Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

On 28 Jun 2021, at 23:44, Anup Pandey ***@***.***> wrote: Is anyone aware of following error while using GAP with LAMMPS? I built LAMMPS with user-quip package following the instructions in a README file (/lib/quip/README) and there was no issue with the built. I tested the exec in the LAMMPS quip tutorial. SYSTEM ABORT: quip_lammps_wrapper: quip_potential not initialised Thank you. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[ap439111]

Here is the full output :
LAMMPS (29 Oct 2020)
using 1 OpenMP thread(s) per MPI task
Reading data file ...
orthogonal box = (0.0000000 0.0000000 0.0000000) to (5.4310000 5.4310000 5.4310000)
1 by 1 by 1 MPI processor grid
reading atoms ...
8 atoms
read_data CPU = 0.001 seconds
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 4.2258
ghost atom cutoff = 4.2258
binsize = 2.1129, bins = 3 3 3
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair quip, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Setting up Verlet run ...
Unit style : metal
Current step : 0
Time step : 0.001
SYSTEM ABORT: quip_lammps_wrapper: quip_potential not initialised

[ap439111]

Hi all,

I have a problem computing elastic constants using GAP with lammps. The model fits well for both the energy and forces (for both training and validation sets). The training data comprise of structures (random multi-component alloys) at different volumes//shapes and also some MD trajectories at different temperatures. The EC calculated using GAP and Lammps are way off the DFT/experimental values. They differ a lot for two random alloys with same compositions (DFT values are close). I was wondering if I have missed something.

Thank you.
Anup

[gabor1]

How good are the stress or virial predictions on your train and test sets? — Gábor Gabor Csanyi Professor of Molecular Modelling Engineering Laboratory University of Cambridge

…

On Wed, Nov 24, 2021 at 18:54, Anup Pandey ***@***.***> wrote: Hi all, I have a problem computing elastic constants using GAP with lammps. The model fits well for both the energy and forces (for both training and validation sets). The training data comprise of structures (random multi-component alloys) at different volumes//shapes and also some MD trajectories at different temperatures. The EC calculated using GAP and Lammps are way off the DFT/experimental values. They differ a lot for two random alloys with same compositions (DFT values are close). I was wondering if I have missed something. Thank you. Anup — You are receiving this because you commented. Reply to this email directly, [view it on GitHub](#305 (reply in thread)), or [unsubscribe](https://github.com/notifications/unsubscribe-auth/ABIANZYGIZSW5GHRW46MX3DUNUYERANCNFSM457GGT6Q). Triage notifications on the go with GitHub Mobile for [iOS](https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675) or [Android](https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub).

[ap439111]

I would say reasonable. I am adding some more data to improve the fit though. Btw I have used the viral (total in the attached figure) from the vasp OUTCAR. It is in the correct format for GAP?

Thank you.

Anup

[ap439111]

Virials for train and validation sets

:

[mcaroba]

Are you parsing VASP stress tensor in any way? VASP reports it with the wrong sign but there is a minus the volume factor going from the stress tensor to the virial. To check for consistency, you can try this with ASE:

from ase.io import read,write
atoms = read("OUTCAR")
atoms.info["virial"] = -atoms.get_volume() * atoms.get_stress(voigt=False)
write("atoms.xyz", atoms)

Does the virial array in the atoms.xyz file look like you'd expect? If not, there could be an inconsistency on how you're feeding virial data to gap_fit.

[ap439111]

I have just tested your script and the virials I have used in my fitting matches the one from your script. So, it is actually the one I pointed out above (the line above "in kb"). It was one of my biggest concerns, thank you for the script.

[gabor1]

you have significant overfitting in your virial data. (what about your forces and energies, are they overfitted too?). are you plotting the total virial or the virial per atom? basic questions: have you triple checked your units and sign? we require the virial in eV, vast provides stresses on kbar (I think). the virial sigma (regularisation parameter) is per atom, so for small cells you might want to use 0.1 or 0.05. For large cells, the virial data is not very useful, you can even skip it. Where you have very large virials, i.e. for highly distorted cells, you might want to increase to virial sigmas quite a bit, lest they compromise your fit of the small virials that enter the elastic constants directly. as a sanity check, you could remove large virials from the training data and see if that gives you better elastic constants. how are you generating the distorted cells? one very good method is to do them essentially randomly on the lattice vectors, using e.g. a high temperature Monte Carlo. a very bad method is to pick a few directions and displace the cell multiple times along fixed directions. — Gábor Gabor Csanyi Professor of Molecular Modelling Engineering Laboratory University of Cambridge ‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐

…

On Wednesday, November 24th, 2021 at 20:08, Anup Pandey ***@***.***> wrote: Virials for train and validation sets [Screen Shot 2021-11-24 at 1 03 27 PM](https://user-images.githubusercontent.com/18051735/143306539-22297ba0-b5fd-4340-9c45-95542d47fbbb.png) : [Screen Shot 2021-11-24 at 1 03 45 PM](https://user-images.githubusercontent.com/18051735/143306545-11380866-e6fa-43e8-a771-b96f86c7ffe0.png) — You are receiving this because you commented. Reply to this email directly, [view it on GitHub](#305 (reply in thread)), or [unsubscribe](https://github.com/notifications/unsubscribe-auth/ABIANZZZIFUIJPZTPOBNQQ3UNVA2HANCNFSM457GGT6Q). Triage notifications on the go with GitHub Mobile for [iOS](https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675) or [Android](https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub).

[ap439111]

I am adding more data to correct the over fitting. When I trained the model with only energy and forces, it does really well with both the energy and forces (also with EV plot compared to the DFT results). Then I added virials to calculate the elastic constants (I am not sure if I need the virials for that, isn't energy and forces sufficient?). Is there a rule of thumb to consider a cell as "large" and skip the virials?

I have randomly distorted the cell with the maximum distortion in angles not to exceed 2 degrees from orthogonal. Then vasp is used for energy/forces calculations.

Thank you for your suggestions.

[gabor1]

one very lazy way to test the virial data is to not fit on virials at all, and then plot GAP virials against DFT virials. they won't be great, but any unit conversion problem or sign problem will be immediately apparent. — Gábor Gabor Csanyi Professor of Molecular Modelling Engineering Laboratory University of Cambridge ‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐

…

On Wednesday, November 24th, 2021 at 23:02, Anup Pandey ***@***.***> wrote: I have just tested your script and the virials I have used in my fitting matches the one from your script. So, it is actually the one I pointed out above (the line above "in kb"). It was one of my biggest concerns, thank you for the script. — You are receiving this because you commented. Reply to this email directly, [view it on GitHub](#305 (reply in thread)), or [unsubscribe](https://github.com/notifications/unsubscribe-auth/ABIANZ5LLOCUDRSRTQ42QJTUNVVGVANCNFSM457GGT6Q). Triage notifications on the go with GitHub Mobile for [iOS](https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675) or [Android](https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub).

[gabor1]

have you checked how good the virials were before you added virial data? and how good the elastic constants were? if an EV curve is good, the curvature of that is basically the bulk modulus which is a combination of elastic constants. virials (correctly regularised) do help elastic constants. I would consider a cell "small" from this perspective it has < 10 atoms. of course if you have inconsistent virial data, then that makes the fit worse. how is your k-point sampling? stresses are notoriously hard to converge, we typically use a k-point spacing of 0.2 inverse Angstroms or less with VASP. — Gábor Gabor Csanyi Professor of Molecular Modelling Engineering Laboratory University of Cambridge ‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐

…

On Wednesday, November 24th, 2021 at 23:11, Anup Pandey ***@***.***> wrote: I am adding more data to correct the over fitting. When I trained the model with only energy and forces, it does really well with both the energy and forces (also with EV plot compared to the DFT results). Then I added virials to calculate the elastic constants (I am not sure if I need the virials for that, isn't energy and forces sufficient?). Is there a rule of thumb to consider a cell as "large" and skip the virials? I have randomly distorted the cell with the maximum distortion in angles not to exceed 2 degrees from orthogonal. Then vasp is used for energy/forces calculations. Thank you for your suggestions. — You are receiving this because you commented. Reply to this email directly, [view it on GitHub](#305 (reply in thread)), or [unsubscribe](https://github.com/notifications/unsubscribe-auth/ABIANZ6KRWYSUZPBBDPSHHTUNVWJBANCNFSM457GGT6Q). Triage notifications on the go with GitHub Mobile for [iOS](https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675) or [Android](https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub).

[ap439111]

I have 32 atoms cell, so I think I can skip the virials.

I have used the virials from vasp to fit the GAP. Its too expensive to calculate elastic constants for the entire training sets but for the few structures that I have calculated, the EC seems to match well with the other theoretical calculations. Same parameters (tight k-sampling the cut-off energy) are used to generate the training data. Btw if the model is trained from virials, I guess the lammps uses the predictions from the GAP instead from calculating from the forces ?

[gabor1]

well, if all your cells are 32 atoms, then it'll be hard to fit the virials from that, for sure, though not impossible. I still think you are overfitting, what was your virial sigma? maybe too small? virials and forces are very different things. the force is the derivative of the total energy with respect to the atomic coordinates, the virial is the derivative with respect to the deformation tensor. imagine a perfectly symmetrical unit cell: all the forces are exactly zero, but the virials are not unless the lattice vectors have been optimised. I think what you re referring to when you say "calculate from forces" is the special case of a pair potential, for which you can obtain the virial stress from the the pairwise forces. but that isn't true for a many-body potential. — Gábor Gabor Csanyi Professor of Molecular Modelling Engineering Laboratory University of Cambridge ‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐

…

On Wednesday, November 24th, 2021 at 23:27, Anup Pandey ***@***.***> wrote: I have 32 atoms cell, so I think I can skip the virials. I have used the virials from vasp to fit the GAP. Its too expensive to calculate elastic constants for the entire training sets but for the few structures that I have calculated, the EC seems to match well with the other theoretical calculations. Same parameters (tight k-sampling the cut-off energy) are used to generate the training data. Btw if the model is trained from virials, I guess the lammps uses the predictions from the GAP instead from calculating from the forces ? [image](https://user-images.githubusercontent.com/18051735/143324746-cc414238-f8bd-491b-ac64-8c6162fd1788.png) — You are receiving this because you commented. Reply to this email directly, [view it on GitHub](#305 (reply in thread)), or [unsubscribe](https://github.com/notifications/unsubscribe-auth/ABIANZZYV522YILA75KUFRDUNVYFFANCNFSM457GGT6Q). Triage notifications on the go with GitHub Mobile for [iOS](https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675) or [Android](https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub).

[ap439111]

So I need to include virials for ECs. The sigma I have used is 0.02.

[gabor1]

yeah I think that is pretty optimistic. try a virial sigma of 0.05 or even 0.1 — Gábor Gabor Csanyi Professor of Molecular Modelling Engineering Laboratory University of Cambridge ‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐

…

On Wednesday, November 24th, 2021 at 23:44, Anup Pandey ***@***.***> wrote: So I need to include virials for ECs. The sigma I have used is 0.02. — You are receiving this because you commented. Reply to this email directly, [view it on GitHub](#305 (reply in thread)), or [unsubscribe](https://github.com/notifications/unsubscribe-auth/ABIANZ3NPV4D7RYC4DAYAV3UNV2HNANCNFSM457GGT6Q). Triage notifications on the go with GitHub Mobile for [iOS](https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675) or [Android](https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub).

[ap439111]

I will do that. Thank you.

[ap439111]

with viral sigma 0.05, I think there is no overfitting. But still most of the EC components are way off. I have trained on ~1000 structures. Is it the issue of data ?

s.

[jameskermode]

Check carefully the ordering of the elements in the virial tensor: Vasp outputs them in an unusual order which does not follow the Voigt convention. ASE parses this correctly, but if you are using your own scripts that could be the source of the discrepancy.

[ap439111]

The order should be good, I have used ASE. I have one confusion, should the single atom energy go at the end of the train.xyz or it can be provided as e0={Mo:-0.392:Nb:-0.750:Ta:-2.374:W:-1.862

}

[gabor1]

Looking at this, yes, it’s very suspicious that the ordering is wrong. your xx and zz are very well correlated. remind me what this system is? why are your YY stress values so much smaller? also, your xx and zz values are very high, asking for 0.05 accuracy for something that is size 50 is 0.1% accuracy, quite big demand! and you don’t get it,you only have 2-3% accuracy even in xx and zz. so I would try to raise the virial sigma to 0.1 as well. have you tried predicting these virials WITHOUT virial data? if you get better predictions than with virial data, then that points to a data problem!

…

-- Gábor Gábor Csányi Professor of Molecular Modelling Engineering Laboratory, University of Cambridge Pembroke College Cambridge Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

On 1 Dec 2021, at 17:12, James Kermode ***@***.***> wrote: Check carefully the ordering of the elements in the viral tensor: Vasp outputs them in an unusual order which does not follow the Voigt convention. ASE parses this correctly, but if you are using your own scripts that could be the source of the discrepancy. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android.

[ap439111]

Thank you for pointing out YY, I was plotting the wrong components. Now it looks good for XX, YY,ZZ.

[gabor1]

you can provide e0, but UNLESS you have the single atoms as part of the train.xyz, you will not get accurate predictions for single atoms.

…

-- Gábor Gábor Csányi Professor of Molecular Modelling Engineering Laboratory, University of Cambridge Pembroke College Cambridge Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

On 1 Dec 2021, at 17:34, Anup Pandey ***@***.***> wrote: The order should be good, I have used ASE. I have one confusion, should the single atom energy go at the end of the train.xyz or it can be provided as e0={Mo:-0.392:Nb:-0.750:Ta:-2.374:W:-1.862 } — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android.

[gabor1]

ok, so what’s happenning now is the the accuracies are quite uniform, around 1-2 eV/atom, and this is quite good for xx-yy-zz, 2%, but very bad fractionally for the off-diagonals because those are smaller, so it’s about 30% accurate only. what happenned from last week? you had about half the error, even for off-diagonals. but let’s think about the ECs. you say they are way off. if I distort the relaxed cell without relaxing atoms, the value of the stress that arises *is* roughly EC * distortion, so there is a contradiction here. presumably when you say ECs are bad, you are measuring them by distorting the cell and relaxing the atoms? If this is the case, are you TRAINING with such configurations? i.e. all those virials that you are showing here (and pretty good accuracy at least for the diagonal parts) are they configurations with distorted unit cells and unrelaxed, relaxed or perturbed atoms? you won’t find good relaxation of atoms in complex unit cells if you train with atoms in their original fractional position.

…

-- Gábor Gábor Csányi Professor of Molecular Modelling Engineering Laboratory, University of Cambridge Pembroke College Cambridge Pembroke College supports CARA. A Lifeline to Academics at Risk. http://www.cara.ngo/

On 1 Dec 2021, at 18:03, Anup Pandey ***@***.***> wrote: Thank you for pointing out YY, I was plotting the wrong components. Now it looks good for XX, YY,ZZ. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android.

[ap439111]

That makes lot more sense. I have distorted the relaxed cell and then measured the ECs without atomic relaxations (or with the atoms at the same fractional positions). Thank you very much.