Using the binding energy to fit the LJ parameter

GULP
1

hello,everyone!
I am new to GULP, I have tried to read the help documentation and some similar topics, but I still have some problems that are not solved, hope someone can help me, thanks.
I intend to use binding energy to fit the non-bonding LJ parameter values of CO2 and HIS (an amino acid). Here’s my code:

fit molecule conv c6 noflag
title
fit the LJ parameter of amine and CO2
end

CO2

cartesian angs
C1 core 5.63 0.0 5.63 0.7
O1 core 6.956118 0.0 5.630001 -0.35
O1 core 4.303882 0.0 5.629999 -0.35
observable
energy ev
-198.8902186 100.0
end

HIS

cartesian angs
C2 core 2.876845 3.391763 6.747024 0.5973
C3 core 4.256139 2.923692 3.234701 0.1868
C4 core 3.477251 2.840953 2.112461 -0.2207
C5 core 3.487032 3.864016 5.450355 -0.0581
C6 core 5.480809 3.601588 1.601743 0.1635
C7 core 3.855906 2.592599 4.627673 -0.0074
O2 core 1.712044 3.072720 6.869164 -0.5679
O2 core 3.768736 3.292290 7.747725 -0.5679
N1 core 4.598141 4.791716 5.665431 -0.4157
N2 core 4.273023 3.277191 1.080208 -0.2795
N3 core 5.508380 3.405782 2.904857 -0.5432
H1 core 5.000993 4.996860 4.755388 0.2719
H2 core 5.335174 4.307586 6.170423 0.2719
H3 core 4.013745 3.331353 0.108037 0.3339
H4 core 2.464974 2.514912 1.959017 0.1862
H5 core 2.696530 4.385351 4.912689 0.1360
H6 core 6.291875 3.966570 0.994294 0.1435
H7 core 4.670754 2.080512 5.145256 0.0367
H8 core 2.992393 1.925867 4.604303 0.0367
H9 core 3.316623 2.927897 8.528148 0.296
observable
energy ev
-580.3261568 100.0
end

the adsorption structure 2

cartesian angs
C1 core 4.462295 6.069737 6.457697 0.7
O1 core 3.279936 6.336698 6.537693 -0.35
O1 core 5.545230 6.527773 6.764285 -0.35
C2 core 2.860954 3.298195 6.766230 0.5973
C3 core 4.189557 2.925964 3.202874 0.1868
C4 core 3.470645 2.780994 2.048981 -0.2207
C5 core 3.476066 3.769916 5.468797 -0.0581
C6 core 5.432760 3.689669 1.621090 0.1635
C7 core 3.773438 2.546866 4.578793 -0.0074
O2 core 1.719451 2.902884 6.832240 -0.5679
O2 core 3.721648 3.296266 7.790850 -0.5679
N1 core 4.699440 4.571213 5.690981 -0.4157
N2 core 4.275929 3.273304 1.052010 -0.2795
N3 core 5.416328 3.500974 2.927335 -0.5432
H1 core 5.155525 4.744034 4.788464 0.2719
H2 core 5.385093 4.059420 6.247361 0.2719
H3 core 4.053827 3.308988 0.068322 0.3339
H4 core 2.493315 2.374708 1.859877 0.1862
H5 core 2.751200 4.414809 4.981404 0.1360
H6 core 6.239465 4.115624 1.047545 0.1435
H7 core 4.549445 1.950342 5.065232 0.0367
H8 core 2.864433 1.949392 4.530873 0.0367
H9 core 3.277900 2.940241 8.580134 0.296
observable
energy ev
-779.5360063 100.0
end

the adsorption structure 2

cartesian angs
C1 core 5.040675 5.980178 6.440142 0.7
O1 core 3.923994 6.455064 6.423932 -0.35
O1 core 6.161825 6.236971 6.825468 -0.35
C2 core 2.986950 3.556547 6.699808 0.5973
C3 core 4.183205 2.906178 3.136584 0.1868
C4 core 3.415524 2.763310 2.013921 -0.2207
C5 core 3.714789 3.863505 5.409641 -0.0581
C6 core 5.392414 3.589631 1.497416 0.1635
C7 core 3.811396 2.590440 4.540752 -0.0074
O2 core 1.798085 3.338302 6.729635 -0.5679
O2 core 3.789410 3.509201 7.768559 -0.5679
N1 core 5.049302 4.448606 5.648516 -0.4157
N2 core 4.206164 3.198228 0.978146 -0.2795
N3 core 5.415755 3.433988 2.805624 -0.5432
H1 core 5.530160 4.553344 4.748960 0.2719
H2 core 5.630667 3.829533 6.208965 0.2719
H3 core 3.958828 3.206232 0.002115 0.3339
H4 core 2.416897 2.395070 1.860858 0.1862
H5 core 3.133234 4.617241 4.889107 0.1360
H6 core 6.192320 3.973109 0.890596 0.1435
H7 core 4.547012 1.920042 4.993123 0.0367
H8 core 2.841219 2.095702 4.558024 0.0367
H9 core 3.262208 3.284613 8.554192 0.296
observable
energy ev
-779.5320175 100.0
end

binding energy

observables
reaction
3 -0.3196309 10.0 1 -1.0 2 -1.0 3 1.0
end
observables
reaction
3 -0.3156421 20.0 1 -1.0 2 -1.0 4 1.0
end

vary
shift
end

the lj paramter, the all element of HIS and CO2

lennard epsilon zero 12 6 inter
C1 core C2 core 0.069 3.100 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core C3 core 0.117 3.214 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core C4 core 0.069 3.100 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core C5 core 0.078 3.100 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core C6 core 0.069 3.100 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core C7 core 0.078 3.100 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core O2 core 0.108 2.880 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core N1 core 0.097 3.025 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core N2 core 0.097 3.025 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core N3 core 0.097 3.025 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core H1 core 0.030 1.935 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core H2 core 0.030 1.935 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core H3 core 0.030 1.935 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core H4 core 0.029 2.655 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core H5 core 0.030 2.636 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core H6 core 0.029 2.611 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core H7 core 0.030 2.725 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core H8 core 0.030 2.725 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
C1 core H9 core 0.030 2.636 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core C2 core 0.117 3.214 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core C3 core 0.117 3.214 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core C4 core 0.117 3.214 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core C5 core 0.132 3.214 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core C6 core 0.117 3.214 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core C7 core 0.132 3.214 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core O2 core 0.183 2.994 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core N1 core 0.165 3.139 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core N2 core 0.165 3.139 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core N3 core 0.165 3.139 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core H1 core 0.050 2.049 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core H2 core 0.050 2.049 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core H3 core 0.050 2.049 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core H4 core 0.049 2.769 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core H5 core 0.050 2.750 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core H6 core 0.049 2.725 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core H7 core 0.050 2.839 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core H8 core 0.050 2.839 0.000 12.000 1 0
lennard epsilon zero 12 6 inter
O1 core H9 core 0.050 2.750 0.000 12.000 1 0

output

dump every HIS_CO2.grs

This code is written with reference to example22 and a few similar topics, but I’m not sure about a few things:

  1. The structure above is calculated from aperiodic systems, so I did not add the cell part, is this correct?
  2. As for the LJ parameter, I want to fit the non-bonding parameter, so I only write the inter part, my intuition tells me that this is not correct, but I don’t know how to modify it. Do I need to add other potential functions?
  3. The LJ parameter here, I’ve written down the interaction between every element of HIS and every element of CO2, is that correct?
  4. I know I don’t have enough configurations in my code right now, is that relevant?

I tried to run the above code, but in the result, all the values related to energy are 0, which indicates that my code is faulty, but I do not know how to modify it.
I sincerely hope that you can help me solve this problem, thank you

2

Here are some quick answers:

  • The structure above is calculated from aperiodic systems, so I did not add the cell part, is this correct?

As the output shows, yes this is correct.

  • As for the LJ parameter, I want to fit the non-bonding parameter, so I only write the inter part, my intuition tells me that this is not correct, but I don’t know how to modify it. Do I need to add other potential functions?

If fitting the binding energy then only including the intermolecular part would be OK. Of course if you want the forces on the atoms or to relax the structures then you’ll need to intramolecular, so best just to include it, especially if there is some non-zero contribution from the intramolecular part because the geometry isn’t optimised.

  • The LJ parameter here, I’ve written down the interaction between every element of HIS and every element of CO2, is that correct?

This is part of the design of the force field, so what is correct is what was used in the fitting of the force field.

  • I know I don’t have enough configurations in my code right now, is that relevant?

Very relevant - you don’t have nearly enough data to fit so many parameters.

Some other comments:

  1. Don’t use the c6 keyword since you don’t have periodic boundary conditions, so it won’t do anything.
  2. Your binding energies are zero because the molecules are too close and are covalently bonded & so there is no intermolecular part. You need to exclude bonds or change the covalent radii.
  3. There’s no point in adding an energy shift here and fitting it. As you have different molecules you’d need a shift for each, but there is no point fitting the energies of the molecules (especially to the large values you have) since they will be zero with no intramolecular interaction anyway. As you’re looking at binding energy, only the difference between the configurations matters.
1 Like
3

Thank you for answering my questions.Your answer is very helpful to me.
I have a few more questions that I would like your help with.

  1. I understand that there are currently only a few configurations in my code, but I wanted to check if my code is correct, so I ran it and got the following results:
    ------------------ results -------------------------------
    Start of fitting :

Cycle: 0 Sum sqs: ************** Gnorm: 0.000000 CPU: 0.007
** Hessian calculated **
**** Fit completed successfully ****

Normally, the results here should show the outcomes of multiple iterations, but I don’t understand why the results appear this way. Is this directly related to the fact that my configurations are too few, or could there be other reasons?

  1. In your response, you mentioned, ‘As you’re looking at binding energy, only the difference between the configurations matters.’ Does this mean that the energies of CO2 and HIS do not need to be set as observed values? In that case, is it unnecessary to set the energies of CO2, HIS, and the HIS-adsorbed CO2? Should we only set the energy for the reaction part?

  2. In my understanding, CO2 and HIS are obtained from DFT calculations (geometry optimizations) in two separate files (referred to as files A and B). Subsequently, I combine CO2 and HIS into a new file (referred to as file C) and perform a DFT calculation to obtain the structure. My question is: Should the coordinates of CO2 and HIS be taken from the results in files A and B, or should they be the coordinates from file C before its geometry optimization?

  3. For the final question, when adding configurations, should I include the configurations for CO2, HIS, and CO2 + HIS, or is it sufficient to only include the configuration for CO2 + HIS?

By the way, CO2 + HIS refers to the configuration where the two are combined together.

I’m very sorry for taking up so much of your time with my questions. I want to express my gratitude once again and wish you all the best in your work.

4

Here are some quick answers, though I would recommend talking to your supervisor about the design of the project since many of the answers will depend on what you’re actually trying to achieve:

  1. This tells you that the parameters you are fitting are not actually able to fit your observables as the Gnorm is zero (i.e. changing the parameters doesn’t alter the sum of squares).
  2. Normally for a binding energy the reference is the optimised configurations of the individual species. These energies don’t need to be in the fit - the binding energy just has to be expressed relative to them.
  3. This is one where you need to talk to your supervisor since choice of geometries depends on what you’re looking to achieve.
  4. See point 2 - you only need to include the individual molecules explicit if you’re fitting something that would change their energies as well as the binding energy.