Force field fitting for MOFs (ERROR : Incorrect potential coefficient input)

Dear Users,

I am totally new to GULP, and am currently using it to reparameterize the UFF4MOF force field under certain pressure. the MOF is HKUST-1, and 6 atoms in the asymmetric unit.

The observables are the experimental volume, bond and angle. I followed the example 73 to prepare my input file, and I got this error:

!! ERROR : Incorrect potential coefficient input
!! : Error is apparently on line 81
!!!

Program terminated by processor 0 in potword21

all parts of the interatomic potentials (LJ, bond, angle, dihedral and improper) need to be fitted. here is just a test input file to make GULP run properly. below is my input file:
#---------------------------------------------------------------
fit conp molecule

cell
26.3 26.3 26.3 90.0 90.0 90.0

Fractional coordinates

fractional
Cu core 0.71480 0.28520 0.5000
O core 0.68296 0.24335 0.44768
C1 core 0.69985 0.13474 0.36526
H core 0.7287 0.1202 0.3798
C2 core 0.67809 0.17809 0.38673
C3 core 0.70348 0.20348 0.43100
space
225

anisotropic_pressure 0.00010132 0.00010132 0.00010132 0.00010132 &
0.00010132 0.00010132 # unit: GPa 1 atm = 0.00010132 GPa

observables
volume
18199.0 # unit: Angstroms**3
fbond
1 2 1.951 # Cu-O Angstroms
fbond
1 1 2.619 # Cu-Cu
fbond
5 6 1.499 # C2-C3
fbond
2 6 1.258 # C3-O
fbond
3 4 0.932 # C1-H1
fbond
3 5 1.395 # C1-C2
fangle
1 2 6 122.14 # Cu-O-C3
fangle
5 6 2 116.80 #C2-C3-O1
fangle
3 5 6 119.71 # C1-C2-C3
fangle
5 3 4 120.22 # C2-C1-H
fangle
1 1 2 84.52 # Cu-Cu-O
fangle
2 1 2 89.71 # O-Cu-O
fangle
2 6 2 126.39 # O-C3-O
fangle
5 3 5 119.56 # C2-C1-C2
end

epsilon kcal
O 0.060000 3.118146 1 0 # All below parameters need to be fitted, do it one by one
Cu 0.005000 3.113691 0 0
C1 0.105000 3.430851 0 0
C2 0.105000 3.430851 0 0
C3 0.105000 3.430851 0 0
H 0.044000 2.571134 0 0

lennard epsilon geometric kcal all 0.0 12.0

harm inter bond ##(THIS IS LINE 81, WHERE ERROR OCCURS)##

Cu O 221.999373 2.029549 0.0 # O-Cu
O C3 714.557442 1.269010 0.0 # O C3
Cu Cu 70.296579 3.032715 0.0 # Cu Cu
C1 H 357.440381 1.081418 0.0 # C H
C1 C2 462.655054 1.379256 0.0 # C1-C2
C2 C3 391.669513 1.458000 0.0 # C2-C3

three cosine kcal
Cu O C3 77.210179 0 0 # Cu-O-C3
O Cu O 111.434476 0 0 # O-Cu-O
Cu Cu O 52.958227 0 0 # Cu-Cu-O
C2 C1 H 57.289016 0 0 # C2-C1-H
C1 C2 C3 111.297508 0 0 # C1-C2-C3
C2 C1 C2 102.183280 0 0 # C2-C1-C2
O C3 O 206.778426 0 0 # O-C3-O
C2 C3 O 137.714246 0 0 #C2-C3-O1

uff4 bond kcal
Cu O C3 O 6.737110 2 10.0 10.0 10.0 10.0 0
H C1 C2 C3 3.368555 2 10.0 10.0 10.0 10.0 0
C1 C2 C3 O 1.250000 2 10.0 10.0 10.0 10.0 0

uffoop bond kcal
H C1 C2 C3 2.000000 1.000000 -1.000000 10.0 10.0 10.0

print fit
print final
output xyz
output forcefield
#---------------------------------------------------------------
I also tried to add this part in the middle, but got the same error:
Species
Cu core Cu4+2
O core O_2
C1 core C_R
C2 core C_R
C3 core C_R
H core H_
library …/gulp/gulp-6.2/Libraries/uff4mof

Can anyone give me some suggestions about why the error occurred in the line “harm inter bond” and how to fix it?

Thank you in advance.

wang

The error isn’t actually coming from that line (there may be something in your file that messes up the line count) since this error is from potential input (which would be the line after). You can always check what is expected for each potential using the help.txt (or html version) that comes with GULP, including how many fitting flags. If you look at the harmonic potentials after the line you’ve highlighted then you’ll see that there are 2 fitting flags missing on each of the next 6 lines.

Dear Julian,

Thanks for your reply. then the error could have occurred because of the wrong coefficient settings for the following types of interatomic potential.

After correcting the fitting flags you mentioned, the error is still present. I found something wrong with the angle setting part.

I am repeating what I set in lammps to optimize the parameters, in which the angle style is cosine/periodic (E(three) = (2K/n2) * (1 - cos(n*theta)) in lammps), but I checked from the GULP manual, the closest one is equatorial (E(three) = (2K/n2) * (1 - cos(ntheta)) + 2Kexp(-beta*(r13 - r0))). however, the last term 2Kexp(-beta(r13 - r0)) is impossible to be 0, how can I solve this problem?

wang

Hi Wang,
I didn’t mean to imply that fixing the flags on the harmonic would solve everything, but it allows you to move on to the next error, if the same mistake was present later in the file.
For the angle-bending potential you’ve missed the “lin3” potential which has the same form as LAMMPS as you’ve listed.
Regards,
Julian