Usage

To calculate potential energy surface refer to PotentialCalculation.Ones you have potential energy calculated you can open it for fitting by using

using PotentialFitting

# There is an example potential in test/data directory
using PotentialDB
r = defaultregistry()
data=loadpotential(r,"4")

Dict{String, Any} with 6 entries:
  "Method"   => "CCSD(T)-F12/RI"
  "cluster2" => Cluster{AtomOnlySymbol} of 1 atoms
  "Energy"   => [0.0172099 0.0232891 … 0.0210108 0.0160872; 0.0104979 0.0147825…
  "Points"   => Cluster{AtomOnlySymbol}[Cluster{AtomOnlySymbol} of 6 atoms Clus…
  "cluster1" => Cluster{AtomOnlySymbol} of 5 atoms
  "Basis"    => "cc-pVDZ-F12 cc-pVDZ-F12-CABS cc-pVTZ/C TIGHTSCF"

Potential can be viewed by energy

using GLMakie # or WGLMakie

# Collumn 1 from data
plot_potential(data)

Geometry can be plotted too by calling

# Collumn 1 from data
scan_potential(data)

External programs can also be called to visualize geometry (VMD here)

visualize_points(data["Points"]; command="vmd")

Setting up Molecules

Next part in defining topology for the potential. This is started by creating two molecules. The information is in the loaded file.

m1, m2 = get_molecules(data)

(Molecule of 5 atoms, Molecule of 1 atoms)

If needed atoms can be flagged as identical.

# Atoms 2 and 3 are identical
makeidentical!(m1, (2,3))

Identical(Set[Set([1]), Set([2, 3]), Set([4]), Set([5])])

Potential Topology

Next we need to define topology for the potential.

mpp = MoleculePairPotential(m1,m2, LennardJones())

Molecule Pair Potential
Molecule 1 has 5 atoms (5 unique) 
Molecule 2 has 1 atoms (1 unique) 

Potential topology:

LennardJones,  C6=0.0, C12=0.0 :  1C :  6Ar
LennardJones,  C6=0.0, C12=0.0 :  2O :  6Ar
LennardJones,  C6=0.0, C12=0.0 :  3O :  6Ar
LennardJones,  C6=0.0, C12=0.0 :  4H :  6Ar
LennardJones,  C6=0.0, C12=0.0 :  5H :  6Ar

Fine Tuning Potential

Alternatively potential can be tuned completely by adding potentials one by one.

# Array where topology is saved
topo=[]

#We can push potential to to this array one at the time
push!(topo,
      # Molecule 1 has 5 atoms so index 6 is molecule 2, or argon now
      PairPotentialTopology(LennardJones(), 1,6)
     )

If needed we can specify which atoms should be treated as identical, by adding information for it in the topology.

# Atoms 2 and 3 of molecule 1 have same potential to to atom 1 of molecule 2
push!(
      topo,
      PairPotentialTopology(LennardJones(), [(2,6), (3,6)])
)

If default form of potential is not enough it can be tuned, by giving it as an input.

push!(
      topo,
      PairPotentialTopology(GeneralPowers(-6,-12), 4,6)
)
push!(
     topo,
     PairPotentialTopology(GeneralPowers(-6,-8, -10, -12), 5,6)
)

Here we used general polynomial potential GeneralPowers to make customized polynomial potential.

We can now create potential.

mpp1=MoleculePairPotential(m1,m2)
mpp1.topology = topo

show(mpp1)

Molecule Pair Potential
Molecule 1 has 5 atoms (5 unique)
Molecule 2 has 1 atoms (1 unique)

Potential topology:

LennardJones,  C6=0.0, C12=0.0 :  1C :  6Ar
LennardJones,  C6=0.0, C12=0.0 :  2O :  6Ar
LennardJones,  C6=0.0, C12=0.0 :  3O :  6Ar
r^-6  r^-12   *** 0.0  0.0   :  4H :  6Ar
r^-6  r^-8  r^-10  r^-12   *** 0.0  0.0  0.0  0.0   :  5H :  6Ar

Preparing Data for Fitting

To do fitting itself we need to prepare fit data.

fdata = FitData(mpp, data["Points"], data["Energy"])

# Also this works and can be used add data from different sources
# fdata = FitData(mpp, data1, data2,...)

At this point we can add weights to data.

# If energy is more than 1500 cm⁻¹ weight is zero
setweight_e_more!(fdata, 0, 1500u"cm^-1")

# If energy is less than 80 cm⁻¹ weight is 4
setweight_e_less!(fdata, 4, 80u"cm^-1")

Fitting Potential

We also need to create fitting model. At the current moment only linear models can be used. Here we take normal linear regression, but any linear model supported by MLJ that supports weights can be used.

using MLJ
LinearRegressor = @load LinearRegressor pkg=GLM

solver = LinearRegressor(fit_intercept=false)

[ Info: For silent loading, specify `verbosity=0`.
import MLJGLMInterface ✔

To do fitting itself.

fit_potential!(solver, mpp, fdata)

Molecule Pair Potential
Molecule 1 has 5 atoms (5 unique) 
Molecule 2 has 1 atoms (1 unique) 

Potential topology:

LennardJones,  ϵ=25.317214239685264 cm^-1, σ=3.475573979126618 Å :  1C :  6Ar
LennardJones,  ϵ=67.82972457310365 cm^-1, σ=3.1475292416535487 Å :  2O :  6Ar
LennardJones,  ϵ=67.82972457310365 cm^-1, σ=3.1475292416535487 Å :  3O :  6Ar
LennardJones,  ϵ=5.158672955899481 cm^-1, σ=2.80007665413333 Å :  4H :  6Ar
LennardJones,  ϵ=4.102565181821855 cm^-1, σ=3.1855955525862223 Å :  5H :  6Ar

Inspecting Fitted Potential

You can inspect the fit by calculating RMSD.

plot_potential(data, mpp)

Alternatively you can visualize also geometry at the same time.

scan_potential(data, mpp)