search

ACEsuit / ACE.jl

Parameterisation of Equivariant Properties of Particle Systems
65 stars 15 forks source link

Real EuclideanVector #38

Closed davkovacs closed 3 years ago

davkovacs commented 3 years ago

Currently the ACE.EuclideanVector can only be of type complex, but I think it would make sense to have a real version too.

cortner commented 3 years ago

I agree - this is important, in fact it should be the default. @MatthiasSachs this is your job I'm afraid. Please think about how that ought to be done.

davkovacs commented 3 years ago

Maybe best to comment it here. I set up the basis and the fit the following way: Fit is good if isreal=false fit is very bad if isreal=true in the line pibasis1 = PIBasis(B1p, N, 4.0; property = μ, isreal=false);

al_train = IPFitting.Data.read_xyz(joinpath(@__DIR__, "./water_test.xyz"),
                                  energy_key="", force_key="", dipole_key="dipole")

μ = ACE.EuclideanVector(Complex{Float64}); # site Dipole is an equivariant property

r0 = 1.1;
N = 3;

#degrees
Deg = ACE.TotalDegree(
      Dict( :n => 1.0, :l => 1.5 ),
      Dict( 1 => 1, 2 => 12/10, 3 => 12/8 )
      );

B1p = ACEatoms.ZμRnYlm_1pbasis(species = [:O, :H],
                               r0 = r0, rin = 0.5, rcut = 4.0, pin = 2,
                               D = Deg, maxdeg = 12);

pibasis1 = PIBasis(B1p, N, 4.0; property = μ, isreal=false);
ace_O = ACE.LinearACEModel(SymmetricBasis(pibasis1, μ));
ace_H = ACE.LinearACEModel(SymmetricBasis(pibasis1, μ));

Vsite = ACEatoms.ACESitePotential(Dict(:O => ace_O, :H => ace_H));

Bsite = ACEatoms.basis(Vsite);
@show length(Bsite);

B = Bsite
function JuLIP.MLIPs.combine(_B::typeof(B), coeffs::AbstractVector)
    @assert _B === B
    ACE.set_params!(Vsite, coeffs)
    return Vsite
end

dbname = ""
dB = LsqDB(dbname, Bsite, al_train);

weights = Dict("default" => Dict("MU" => 1.0 ));

IP_dip, lsqinfo = IPFitting.Lsq.lsqfit(dB, 
                                solver=(:lsqr, (0.001, 1e-7)),
                                asmerrs=true, weights=weights);
davkovacs commented 3 years ago

and here is the water_test.xyz

3
Lattice="0.0 0.0 50.0 0.0 50.0 0.0 50.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 config_type=monomer energy=-2050.069116897478 dipole="2.0095555787924324 0.93706566395159 0.0" pbc="F F F"
O       25.00000000      25.00000000      25.00000000    -124.32501084     -57.97362922      -0.00000000
H       25.60000000      25.00000000      25.00000000      89.70712689     -30.08282096       0.00000000
H       25.38567257      25.45962667      25.00000000      34.61788395      88.05645018       0.00000000
3
Lattice="0.0 0.0 50.0 0.0 50.0 0.0 50.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 config_type=monomer energy=-2051.2076537559064 dipole="1.948248654403379 0.9926790762912456 -0.0" pbc="F F F"
O       25.00000000      25.00000000      25.00000000    -120.38189842     -61.33763418       0.00000000
H       25.60000000      25.00000000      25.00000000      88.31849148     -24.53455302      -0.00000000
H       25.35267115      25.48541020      25.00000000      32.06340694      85.87218720       0.00000000
3
Lattice="0.0 0.0 50.0 0.0 50.0 0.0 50.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 config_type=monomer energy=-2052.142505647024 dipole="1.8838916342735765 1.044276529238915 -0.0" pbc="F F F"
O       25.00000000      25.00000000      25.00000000    -116.23218315     -64.42858598       0.00000000
H       25.60000000      25.00000000      25.00000000      87.21445333     -20.28057076      -0.00000000
H       25.31795156      25.50882886      25.00000000      29.01772982      84.70915674      -0.00000000
3
Lattice="0.0 0.0 50.0 0.0 50.0 0.0 50.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 config_type=monomer energy=-2052.919352076845 dipole="1.8171199550156816 1.091832605330092 -0.0" pbc="F F F"
O       25.00000000      25.00000000      25.00000000    -111.90639789     -67.24013068       0.00000000
H       25.60000000      25.00000000      25.00000000      86.33186761     -16.93856440       0.00000000
H       25.28168294      25.52976856      25.00000000      25.57453029      84.17869508       0.00000000
3
Lattice="0.0 0.0 50.0 0.0 50.0 0.0 50.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 config_type=monomer energy=-2053.570640885456 dipole="1.748467383384324 1.135474391887308 -0.0" pbc="F F F"
O       25.00000000      25.00000000      25.00000000    -107.42893711     -69.76510622       0.00000000
H       25.60000000      25.00000000      25.00000000      85.62199800     -14.25068836      -0.00000000
H       25.24404199      25.54812727      25.00000000      21.80693911      84.01579457       0.00000000
3
Lattice="0.0 0.0 50.0 0.0 50.0 0.0 50.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 config_type=monomer energy=-2054.119877447974 dipole="1.6784168512051238 1.1752527238395754 0.0" pbc="F F F"
O       25.00000000      25.00000000      25.00000000    -102.81934761     -71.99487995       0.00000000
H       25.60000000      25.00000000      25.00000000      85.04719808     -12.04193666       0.00000000
H       25.20521209      25.56381557      25.00000000      17.77214953      84.03681662       0.00000000
3
Lattice="0.0 0.0 50.0 0.0 50.0 0.0 50.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 config_type=monomer energy=-2054.584461734382 dipole="1.6073750379101799 1.2112438535804126 0.0" pbc="F F F"
O       25.00000000      25.00000000      25.00000000     -98.09370460     -73.91892842       0.00000000
H       25.60000000      25.00000000      25.00000000      84.57826515     -10.19233059      -0.00000000
H       25.16538241      25.57675702      25.00000000      13.51543945      84.11125901       0.00000000
3
Lattice="0.0 0.0 50.0 0.0 50.0 0.0 50.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 config_type=monomer energy=-2054.9775831919774 dipole="1.5355961181445552 1.2434986160388326 0.0" pbc="F F F"
O       25.00000000      25.00000000      25.00000000     -93.26705316     -75.52611935       0.00000000
H       25.60000000      25.00000000      25.00000000      84.19237442      -8.61826301      -0.00000000
H       25.12474701      25.58688856      25.00000000       9.07467874      84.14438236       0.00000000
3
Lattice="0.0 0.0 50.0 0.0 50.0 0.0 50.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 config_type=monomer energy=-2055.3094762030473 dipole="1.4632071792307806 1.271940758821316 0.0" pbc="F F F"
O       25.00000000      25.00000000      25.00000000     -88.35454883     -76.80542799       0.00000000
H       25.60000000      25.00000000      25.00000000      83.87182394      -7.26062737      -0.00000000
H       25.08350386      25.59416084      25.00000000       4.48272489      84.06605536      -0.00000000
3
Lattice="0.0 0.0 50.0 0.0 50.0 0.0 50.0 0.0 0.0" Properties=species:S:1:pos:R:3:forces:R:3 config_type=monomer energy=-2055.588272513455 dipole="1.3902590560978685 1.2964431946053325 -0.0" pbc="F F F"
O       25.00000000      25.00000000      25.00000000     -83.37250117     -77.74612758       0.00000000
H       25.60000000      25.00000000      25.00000000      83.60291784      -6.07706442       0.00000000
H       25.04185388      25.59853843      25.00000000      -0.23041668      83.82319200      -0.00000000
davkovacs commented 3 years ago

@MatthiasSachs Please let me know if you need any more information, I am afraid I have not managed to observe much more than this so far.

cortner commented 3 years ago

@MatthiasSachs Can you please start by reproducing David's results and then we need to debug this. I'm really surprised, I thought our tests were quite comprehensive.

MatthiasSachs commented 3 years ago

I am working on this. However, al_train = IPFitting.Data.read_xyz(joinpath(@__DIR__, "./water_test.xyz"), energy_key="", force_key="", dipole_key="dipole") throws an error. I think the error is caused by the keyword dipole_key="dipole" I am currently working on the dev branch of IPFitting.

MatthiasSachs commented 3 years ago

@cortner @davkovacs I guess David has implemented that extension on his branch. Is there a way to quickly merge the relevant change into the developer branch of IPFitting?

davkovacs commented 3 years ago

Can you just checkout the dipole branch of IPFitting?

MatthiasSachs commented 3 years ago

@davkovacs could you send me the status of your julia environment / package manager ? I checked out the relevant dipole branches both for ACEatoms and IPFitting, but I am still having trouble getting things running...

davkovacs commented 3 years ago

Sorry about that, I have probably changed JuLIP as well to get this to work: Project.toml:

[deps]
ACE = "3e8ccfd2-c8b0-11ea-32f1-f3a5990fd77a"
ACEatoms = "1e34e032-0b37-4839-a012-196b35408c3c"
ACEbase = "14bae519-eb20-449c-a949-9c58ed33163e"
IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
IPFitting = "3002bd4c-79e4-52ce-b924-91256dde4e52"
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
JuLIP = "945c410c-986d-556a-acb1-167a618e0462"
NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
MatthiasSachs commented 3 years ago

oh, I should have been more precise: could you send me the output of ] status and the branches of all packages in development mode?

davkovacs commented 3 years ago

status:

      Status `~/.julia/environments/v1.6/Project.toml`
  [3e8ccfd2] ACE v0.10.4 `~/.julia/dev/ACE`
  [1e34e032] ACEatoms v0.0.5-pre `~/.julia/dev/ACEatoms`
  [14bae519] ACEbase v0.1.6 `~/.julia/dev/ACEbase`
  [7073ff75] IJulia v1.23.2
  [3002bd4c] IPFitting v0.5.0 `~/.julia/dev/IPFitting`
  [42fd0dbc] IterativeSolvers v0.9.1
  [945c410c] JuLIP v0.12.3 `~/.julia/dev/JuLIP`
  [d9ec5142] NamedTupleTools v0.13.7
  [90137ffa] StaticArrays v1.2.2
  [e88e6eb3] Zygote v0.6.11
cortner commented 3 years ago

that's getting more difficult to reproduce if you've changed JuLIP - can you please explain what you changed there?

davkovacs commented 3 years ago

branches: ACE: 'origin/dev-v0.10.x'. ACEatoms: 'origin/dipole'. IPFitting: 'origin/dipole' JuLIP: 'origin/master'

davkovacs commented 3 years ago

I was wrong there, sorry nothing changed in JuLIP, I was planning to, but you warned me against it.

cortner commented 3 years ago

@MatthiasSachs Try the following configuration for the Project.toml:


[deps]
ACE = "3e8ccfd2-c8b0-11ea-32f1-f3a5990fd77a"
ACEatoms = "1e34e032-0b37-4839-a012-196b35408c3c"
ACEbase = "14bae519-eb20-449c-a949-9c58ed33163e"
IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
IPFitting = "3002bd4c-79e4-52ce-b924-91256dde4e52"
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
JuLIP = "945c410c-986d-556a-acb1-167a618e0462"
NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[compat]
ACE = "0.10"
ACEbase = "0.1.5"
julia = "1.5,1.6"
cortner commented 3 years ago

and then dev IPFitting and ACEatoms both on the dipole branches

My status:

(david) pkg> status
      Status `~/gits/temp/david/Project.toml`
  [3e8ccfd2] ACE v0.10.4
  [1e34e032] ACEatoms v0.0.4 `~/.julia/dev/ACEatoms`
  [14bae519] ACEbase v0.1.6
  [7073ff75] IJulia v1.23.2
  [3002bd4c] IPFitting v0.5.0 `~/.julia/dev/IPFitting`
  [42fd0dbc] IterativeSolvers v0.9.1
  [945c410c] JuLIP v0.12.3
  [d9ec5142] NamedTupleTools v0.13.7
  [90137ffa] StaticArrays v1.2.2
  [e88e6eb3] Zygote v0.6.12
davkovacs commented 3 years ago

I think we have solved this @MatthiasSachs I was putting the isreal=true to the PIbasis and not the SymmetricBasis, sorry about that. All works now the same irrespective of putting real or complex there.

MatthiasSachs commented 3 years ago

Ok, great! @cortner Would it make sense to output a warning if an inconsistent combination of the isreal flags is used in SymmetricBasis and PIBasis?

cortner commented 3 years ago

I'd like to come back to this with the hindsight of a bit more experience. This is an interesting design issue actually.

If you give me a property are you able to tell me whether there should be a real call between the product basis and the symmetrisation, and another real call after symmetrisation? Or would the user have to specify that manually?

IOW, I'd like to make this fully automated, allowing the property to auto-magically insert the call to real before and/or after the symmetrisation is applied. The main design question is whether the product basis needs to know about it, or not. This would save a little bit of performance especially at low correlation order but I'm tempted to just ignore that.

@MatthiasSachs Do you agree this would be the right way forward, then maybe we can work on a PR together?

cortner commented 3 years ago

As I'm thinking about it, I think I want to suggest the following: 

pre_a2b(aa, property)      # apply this before the A2B mapping 
post_a2b(aa, property)    # apply this after the A2B mapping

maybe there are better names...

cortner commented 3 years ago

This should be fully resolved with #51