update snap

src/MDP/Potential/Potential.jl

@@ -11,9 +11,10 @@ module Potential
 
 using Revise, LinearAlgebra
 
-export empiricalpotential, addpotential, getpotential, accumarray, boundingbox, domainmaps
-export atomlist, neighborlist, fullneighborlist, neighsingles, neighpairlist, neighpairs 
-export neightripletlist, neightriplets, neighquadletlist, neighquadlets 
+export empiricalpotential, empiricaldescriptors, lammpspotential, lammpssnapdescriptors
+export addpotential, getpotential, initsna, snappotential, snapdescriptors
+export boundingbox, domainmaps, atomlist, neighborlist, fullneighborlist, neighsingles, accumarray 
+export neighpairlist, neighpairs, neightripletlist, neightriplets, neighquadletlist, neighquadlets 
 export tallysingle, tallypair, tallytriplet, tallyquadlet, findatomtype, periodicimages
 
 include("accumarray.jl");
@@ -34,7 +35,12 @@ include("neightriplets.jl");
 include("neighquadletlist.jl");
 include("neighquadlets.jl");
 include("empiricalpotential.jl");
+include("empiricaldescriptors.jl");
 include("tally.jl");
+include("snastruct.jl");
+include("snapcpp.jl");
+include("lammpspotential.jl");
+include("lammpssnapdescriptors.jl");
 
 end
 

src/MDP/Potential/empiricaldescriptors.jl

@@ -0,0 +1,127 @@
+function empiricaldescriptors(x, q, t, a, b, c, pbc, rcutmax, eta, kappa, potential)
+
+    # compute full neighbor list
+    y, alist, neighlist, neighnum = fullneighborlist(x, a, b, c, pbc, rcutmax);
+    if !isempty(q)
+        q = q[:,alist];
+    end
+
+    dim, N = size(x);
+    d = [];
+    dd = [];    
+    for i = 1:length(potential)        
+
+        style, bondtype, bondedatomtypes, potentialfunc, mu, rcut = getpotential(potential[i]);    
+        potentialfunc = getfield(Main, Symbol(potentialfunc));
+
+        # single potentials
+        if (style == "single") & (bondtype == "nonbonded")            
+            ilist = Array(1:N);                        
+         elseif (style == "single") & (bondtype == "bonded")     
+            typei = bondedatomtypes[1];
+            ilist = findatomtype(Array(1:N), t, typei);      
+        end
+        if (style == "single")
+            xi, qi, ai, ti = neighsingles(x, q, t, ilist);        
+            ei, fi = potentialfunc(xi, qi, ti, mu, eta, kappa);
+            ea, fa = tallysingle(ei, fi, ai);   
+            d = [d; sum(ea)];         
+            if length(dd)==0
+                dd = fa;
+            else
+                dd = cat(dd, fa, dims=3)
+            end    
+        end
+
+        # pair potentials
+        if (style == "pair") & (bondtype == "nonbonded")            
+            ilist = Int64.(Array(1:N));                        
+            pairlist, pairnum = neighpairlist(y, ilist, neighlist, neighnum, rcut*rcut);
+         elseif (style == "pair") & (bondtype == "bonded")     
+            typei = bondedatomtypes[1];
+            typej = bondedatomtypes[2];            
+
+            ilist = findatomtype(Array(1:N), t, typei);      
+            pairlist, pairnum = neighpairlisttypej(y, ilist, alist, t, typej, neighlist, neighnum, rcut*rcut);                                                
+        end
+        if (style == "pair")
+            xij, ai, aj, ti, tj, qi, qj = neighpairs(y, q, pairlist, pairnum, t, ilist, alist);            
+            eij, fij = potentialfunc(xij, qi, qj, ti, tj, mu, eta, kappa);            
+            for m = 1:size(eij,2)
+                ea, fa = tallypair(eij[:,m], fij[:,:,m], ai, aj);            
+                d = [d; sum(ea)];
+                if length(dd)==0
+                    dd = fa;
+                else
+                    dd = cat(dd, fa, dims=3)
+                end    
+            end
+        end
+
+        # triplet potentials
+        if (style == "triplet") & (bondtype == "nonbonded")            
+            ilist = Array(1:N);                     
+            pairlist, pairnum = neighpairlist(y, ilist, neighlist, neighnum, rcut*rcut);        
+         elseif (style == "triplet") & (bondtype == "bonded")     
+            typei = bondedatomtypes[1];
+            typej = bondedatomtypes[2];            
+            typek = bondedatomtypes[3];            
+
+            ilist = findatomtype(Array(1:N), t, typei);      
+            pairlist, pairnum = neighpairlisttypej(y, ilist, alist, t, [typej; typek], neighlist, neighnum, rcut*rcut);                                                
+        end
+        if (style == "triplet")
+            tripletlist, tripletnum = neightripletlist(pairlist, pairnum, ilist);
+            xij, xik, ai, aj, ak, ti, tj, tk, qi, qj, qk = neightriplets(y, q, tripletlist, tripletnum, atomtype, ilist, alist);
+            eij, fij, fik = potentialfunc(xij, xik, qi, qj, qk, ti, tj, tk, mu, eta, kappa);                    
+            for m = 1:size(eij,2)
+                ea, fa = tallytriplet(eij[:,m], fij[:,:,m], fik[:,:,m], ai, aj, ak);                             
+                d = [d; sum(ea)];
+                if length(dd)==0
+                    dd = fa;
+                else
+                    dd = cat(dd, fa, dims=3)
+                end    
+            end
+        end
+
+        # quadlet potentials
+        if (style == "quadlet") & (bondtype == "nonbonded")            
+            ilist = Array(1:N);                        
+            pairlist, pairnum = neighpairlist(y, ilist, neighlist, neighnum, rcut*rcut);        
+         elseif (style == "quadlet") & (bondtype == "bonded")     
+            typei = bondedatomtypes[1];
+            typej = bondedatomtypes[2];            
+            typek = bondedatomtypes[3];      
+            typel = bondedatomtypes[4];      
+
+            ilist = findatomtype(Array(1:N), t, typei);      
+            pairlist, pairnum = neighpairlisttypej(y, ilist, alist, t, [typej typek typel], neighlist, neighnum, rcut*rcut);                                                
+        end
+        if (style == "quadlet")
+            quadletlist, quadletnum = neighquadletlist(pairlist, pairnumsum, ilist);        
+            xij, xik, xil, ai, aj, ak, al, ti, tj, tk, tl, qi, qj, qk, ql = 
+                        neighquadlets(x, q, quadletlist, quadletnum, atomtype, ilist, alist);
+            eij, fij, fik, fil = potentialfunc(xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, mu, eta, kappa);            
+            for m = 1:size(eij,2)
+                ea, fa = tallyquadlet(eij[:,m], fij[:,:,m], fik[:,:,m], fil[:,:,m], ai, aj, ak, al);                            
+                d = [d; sum(ea)];
+                if length(dd)==0
+                    dd = fa;
+                else
+                    dd = cat(dd, fa, dims=3)
+                end    
+            end
+        end
+
+    end    
+
+    return d, dd
+
+end
+
+
+
+
+
+

src/MDP/Potential/lammpspotential.jl

@@ -0,0 +1,81 @@
+using LAMMPS
+
+function lammpspotential(x, t, a, b, c, pbc, unitstyle, pair_style, pair_coeff)
+
+    lmp = LMP()
+
+    # handle boundary conditions
+    bcs = ""
+    for i = 1:3
+        if pbc[i] == 1
+            bcs = bcs * " p"
+        else  
+            bcs = bcs * " f"
+        end
+    end
+
+    command(lmp, "units " * unitstyle)
+    command(lmp, "boundary " * bcs)
+    command(lmp, "atom_style atomic")    
+    command(lmp, "atom_modify map array")
+    command(lmp, "box tilt large")
+
+    # create simulation box
+    xlo = 0.0; ylo = 0.0; zlo = 0.0;
+    xhi = a[1]; yhi = b[2]; zhi = c[3];
+    xy = b[1]; xz = c[1]; yz = c[2];    
+    if (abs(xy) + abs(xz) + abs(yz)) <= 1e-10
+        command(lmp, "region simbox block $xlo $xhi $ylo $yhi $zlo $zhi units box")
+    else
+        command(lmp, "region simbox prism $xlo $xhi $ylo $yhi $zlo $zhi $xy $xz $yz units box")
+    end
+    command(lmp, "create_box 1 simbox")
+
+    # pair potential 
+    command(lmp, "pair_style " * pair_style)
+    for i = 1:length(pair_coeff)
+        command(lmp, "pair_coeff " * pair_coeff[i]) 
+    end
+
+    # set mass 
+    ntypes = length(unique(t))
+    for i = 1:ntypes
+        itype = Int32(i)
+        command(lmp, "mass $itype 1.0")
+    end
+
+    # create atom types and atom positions
+    n = size(x,2)
+    id = Int32.(Array(1:(n)))    
+    v = 0.0*x
+    image = zeros(Int32, n)
+    bexpand = 0
+    LAMMPS.API.lammps_create_atoms(lmp, n, id, Int32.(t), x, v, image, bexpand)
+
+    # run lammps
+    command(lmp, "run 0")
+
+    # potential energy
+    etot = LAMMPS.API.lammps_get_thermo(lmp, "pe")
+    pe = Float64(etot)*n
+
+    # extract forces
+    forces = extract_atom(lmp, "f")
+    id = extract_atom(lmp, "id")
+    pos = extract_atom(lmp, "x")
+    type = extract_atom(lmp, "type")
+
+    return pe, forces, id, pos, type
+
+    # pos = extract_atom(lmp, "x")
+    # t = extract_atom(lmp, "type")
+    # i = extract_atom(lmp, "id")    
+    # nlocal = extract_global(lmp, "nlocal")
+    # ntypes = extract_global(lmp, "ntypes")
+end
+
+
+
+
+
+

src/MDP/Potential/lammpssnapdescriptors.jl

@@ -0,0 +1,123 @@
+using LAMMPS
+
+function lammpssnapdescriptors(x, t, a, b, c, pbc, unitstyle, pair_style, pair_coeff, snapparam, elemradius, elemweight)
+
+    lmp = LMP()
+
+    # handle boundary conditions
+    bcs = ""
+    for i = 1:3
+        if pbc[i] == 1
+            bcs = bcs * " p"
+        else  
+            bcs = bcs * " f"
+        end
+    end
+
+    command(lmp, "units " * unitstyle)
+    command(lmp, "boundary " * bcs)
+    command(lmp, "atom_style atomic")    
+    command(lmp, "atom_modify map array")
+    command(lmp, "box tilt large")
+
+    # create simulation box
+    xlo = 0.0; ylo = 0.0; zlo = 0.0;
+    xhi = a[1]; yhi = b[2]; zhi = c[3];
+    xy = b[1]; xz = c[1]; yz = c[2];    
+    if (abs(xy) + abs(xz) + abs(yz)) <= 1e-10
+        command(lmp, "region simbox block $xlo $xhi $ylo $yhi $zlo $zhi units box")
+    else
+        command(lmp, "region simbox prism $xlo $xhi $ylo $yhi $zlo $zhi $xy $xz $yz units box")
+    end
+    command(lmp, "create_box 1 simbox")
+
+    # pair potential 
+    command(lmp, "pair_style " * pair_style)
+    for i = 1:length(pair_coeff)
+        command(lmp, "pair_coeff " * pair_coeff[i]) 
+    end
+
+    # set mass 
+    ntypes = length(unique(t))
+    for i = 1:ntypes
+        itype = Int32(i)
+        command(lmp, "mass $itype 1.0")
+    end
+
+    # create atom types and atom positions
+    natom = size(x,2)
+    id = Int32.(Array(1:(natom)))    
+    v = 0.0*x
+    image = zeros(Int32, natom)
+    bexpand = 0
+    LAMMPS.API.lammps_create_atoms(lmp, natom, id, Int32.(t), x, v, image, bexpand)
+
+    command(lmp, "compute PE all pe")
+    command(lmp, "compute S all pressure thermo_temp")
+
+    # snap parammeters 
+    ntypes = Int32(snapparam[1]);  
+    twojmax = Int32(snapparam[2]);  
+    rcutfac = snapparam[3];
+    rfac0 = snapparam[4];
+    rmin0 = snapparam[5];
+    bzeroflag = Int32(snapparam[6]);
+    switchflag = Int32(snapparam[7]);
+    quadraticflag = Int32(snapparam[8]);
+    chemflag = Int32(snapparam[9]);
+    bnormflag = Int32(snapparam[10]);
+    wselfallflag = Int32(snapparam[11]);   
+    flags = " bzeroflag $bzeroflag quadraticflag $quadraticflag switchflag $switchflag bnormflag $bnormflag wselfallflag $wselfallflag"
+
+    radius = ""
+    weight = ""
+    nelements = length(elemradius)
+    for i = 1:nelements
+        elemr = elemradius[i]
+        elemw = elemweight[i]
+        radius = radius * " $elemr"
+        weight = weight * " $elemw"
+    end    
+    keys = " rmin0 $rmin0"
+    if (nelements > 1) & (chemflag==1)                
+        keys = keys * " chem $nelements"
+        for i = 1:nelements
+            itype = Int32(i)
+            keys = keys * " $itype"
+        end
+    end    
+
+    # compute snap descriptors
+    command(lmp, "compute snap all snap $rcutfac $rfac0 $twojmax" * radius * weight * flags * keys)    
+
+    # run lammps
+    command(lmp, "thermo_style custom pe")
+    command(lmp, "run 0")
+
+    # all : [num_bispectrum + num_reference] x [num_energy + num_forces + num_virials]    
+    all = extract_compute(lmp, "snap", LAMMPS.API.LMP_STYLE_GLOBAL,
+                                       LAMMPS.API.LMP_TYPE_ARRAY)        
+
+    # extract bispectrum, derivatives, and virials                                       
+    bi = all[1:end-1,1]
+    bd = all[1:end-1,2:(1+3*natom)]
+    bv = all[1:end-1,(2+3*natom):end];
+
+    num_bispectrum = length(bi) 
+    bd = reshape(bd, (num_bispectrum, 3, natom))
+    bd = permutedims(bd, (2,3,1))
+    bv = permutedims(bv, (2,1))
+
+    # extract reference energy, forces, and stresses
+    energy = all[end,1]
+    forces = reshape(all[end,2:(1+3*natom)], (3, natom))    
+    stresses = all[end,(2+3*natom):end];
+
+    return bi, bd, bv, energy, forces, stresses
+end
+
+
+
+
+
+