Added example of Si using force field; updated manual

Author: feifzhou

README.md

@@ -14,7 +14,7 @@ Then:
 Documentation
 ----------------
 
-Check out [manual.pdf](manual.pdf)
+Check out [manual.rst](manual.rst)
 
 
 Get Involved

examples/Si/csld.in-forcefield

@@ -0,0 +1,126 @@
+[structure]
+prim = POSCAR
+sym_tol = 1E-3
+## Either text file containing 3x3 matrix, or 9 numbers for the same matrix
+#epsilon_inf = 2 0 0   0 2 0   0 0 2
+### Either text file containing n_atom matrices, or n_atom*9 numbers for the same matrices
+#born_charge = 1 0 0  0 1 0  0 0 1    -1 0 0  0 -1 0  0 0 -1
+
+#######################################################################
+[model]
+model_type = LD
+cluster_in = clusters.out
+cluster_out = clusters.out
+symC_in = Cmat.mtx
+symC_out = Cmat.mtx
+max_order = 4
+fractional_distance = False
+cluster_diameter = 10.8 6.1 5.0
+## Be VERY careful with these! Do not change unless you understand it. Bond distance cutoff is adjustable
+cluster_filter = lambda cls: ((cls.order_uniq<=2) or (cls.bond_counts(2.9) >=2)) and cls.is_small_in('fcc333/sc.txt')
+
+# polynomial expansion based on bond length
+# must be used together with options "--ldff_step=2"
+[LDFF]
+# number of basis functions
+num_basis = 12
+# list of which orbit of clusters to include in LDFF
+# the nearest-neighbor pair interaction (0=const, 1=point)
+orbit_indices = 3
+# min, max, interval of sampling points
+interpolation_pts = -0.3 0.3 0.01
+# pair basis functions
+basis_2 = lambda l, x: np.polynomial.legendre.legval(x/0.5, np.identity(13)[l+1])
+
+
+#######################################################################
+[training]
+interface = VASP
+corr_type = f
+corr_in = Amat.mtx
+corr_out = Amat.mtx
+fval_in = fval.txt
+fval_out = fval.txt
+traindat1 = fcc333/SPOSCAR fcc333/dir*
+
+
+
+#######################################################################
+[fitting]
+solution_in = solution_all
+solution_out = sol_FF
+nsubset = 5
+holdsize = 0.09
+## 1 FPC 2 FPC sparse 3 split 4 sparse split
+## 5 split+ right preconditioning 6 sparse split + r preconditioning
+## 101 Bayesian CS
+method = 5
+
+# For weight of L1 or L2 regularization
+mulist =  1E-3 1E-4 1E-5 1E-6 1E-7 1E-9
+maxIter = 300
+tolerance = 1E-6
+subsetsize = 0.85
+lambda = 0.5
+uscale_list = 0.03
+
+submodel1 = pairFF 5
+
+[phonon]
+##  -1 (disabled,default), 0 (dipole)
+#nac= 0
+
+## dispersion
+qpoint_fractional = False
+# 'Auto' or something like "[[10,  [0,0,0],'\\Gamma', [0.5,0.5,0.5], 'X', [0.5,0.5,0], 'K']]"
+wavevector = [[25,  [0,0,0],'\Gamma', [0,0.5,0.5], 'X'], [25, [1,0.5,0.5], 'X', [0.75,0.375,0.375], 'K', [0,0,0], '\Gamma', [0.5, 0.5, 0.5], 'L']]
+#wavevector = Auto
+## THz, meV, eV, cm
+unit = THz
+
+## num. of grid points, 3 integers
+dos_grid = 15 15 15
+## Num. of points in DOS
+nE_dos = 400
+## 0 (Gaussian), 1 (Lorentzian), -1 (tetrahedron method)
+ismear= -1
+## width in THz of Gaussian/Lorentzian smearing
+epsilon = 0.05
+pdos = True
+
+
+## debye velocity averaging
+#debye_t_qfrac = 0.03
+#debye_t_v_intgrid = 40 40
+
+## quasi-harmonic thermal properties (requires dos_grid)
+## Temperature in K: begin end increment
+thermal_T_range = 50 800 50
+thermal_out = thermal_out.txt
+
+## export structures to supercells
+#supercell = 2 0 0    0 2 0    0 0 2
+## temperature in K: start, end, increment
+#snapshot_T_range = 100 800 100
+## list of [kp1, kp2, kp3, index_band]
+#modes = 0.125 0.125 0.25 0   0.125 0.125 0.25 1  0.125 0.125 0.25 2
+## amplitude in Angstrom
+#mode_amplitude = 0 0.03
+
+
+[export_potential]
+#save_pot_cell0 = fcc333  -3 3 3   3 -3 3     3 3 -3
+#save_pot_cell1 = prim111 1 0 0   0 1 0    0 0 1
+#combine_improper = false
+## first 3 integers for size of supercell for which Hessian/2nd order FCs are exported, then order of FCs
+export_shengbte = 5 5 5 2 3 4
+
+
+[prediction]
+interface = VASP
+corr_type = f
+corr_in = Amat_pred.mtx
+corr_out = Amat_pred.mtx
+fval_in = fval_pred.txt
+fval_out = fval_pred.txt
+traindat0 = fcc222/POSCAR fcc222/traj*

examples/Si/fit-forcefield.sh

@@ -0,0 +1,13 @@
+SCRIPT=../../scripts/csld_main
+# first fit force field coefficients only
+$SCRIPT --ldff_step 2 -f csld.in-forcefield --phonon_step -1 
+# next read coefficients and fit FCT only
+$SCRIPT --ldff_step 2 -f csld.in-forcefield --phonon_step -1 --symC_step 1 --train_step 1 \
+  --override '[fitting] solution_known = sol_FF' \
+  --override '[fitting] solution_out = sol_LD+FF' \
+  --override '[fitting] submodel1 = FCT 0 1 2 3 4'
+# now the obtained phonon doesn't make sense, but the FF potential should look okay.
+
+# alternatively, 1 fit orders 0,1,2 for phonons; 2. fit FF (order "5"); 3. fit orders 3 4
+# such that the phonon spectrum is correct but the FF looks weird.
+# your choice

manual.rst

@@ -311,6 +311,19 @@ The csld_main script runs in the following steps, each controlled by a command-l
   "solution_in", str , "filename for loading previous solution instead of fitting"
   "solution_known", str, "filename for previously obtained parameter phi_in. The force predicted by phi_in will be subtracted from the total force. This is useful in conjunction with submodel to fit in several steps, e.g. assuming max_order=4, first fit harmonic terms with submodel=harmonic 1 2; solution_out=sol_2nd, then fit anharmonic terms with solution_known=sol_2nd; submodel1=anh 3 4; solution_out=solution_all"
 
+- Pairwise force-field setup step to capture the bulk of the anharmonicity and to make fitting the residual force/energy easier. Turned **off** by default
+
+  - switch "--ldff_step STEP". **0**=off, 2=on.
+
+.. csv-table:: "[LDFF]" section
+  :header: "tag", "value", "description"
+  :widths: 4, 4, 16
+
+  num_basis, int, "Number of basis"
+  orbit_indices, "int [int]", "list of clusters to include in LDFF, usually choose only the nearest-neighbor pairs"
+  interpolation_pts, "int int int", "min, max, interval of sampling points"
+  basis_2, str, "pair basis functions. See Si/csld.in-forcefield for an example"
+
 - Phonon step
 
   - switch "--phonon_step STEP". 0=skip, **1** = compute.