diff --git a/fitpot/makefile.in b/fitpot/makefile.in
index 401be5a..d2b5b58 100644
--- a/fitpot/makefile.in
+++ b/fitpot/makefile.in
@@ -55,7 +55,7 @@ define cd_test
 
 endef
 
-test:
+test: fitpot
 	@$(foreach x,$(wildcard $(fitpot_examples)),$(call cd_test, ${x}))
 
 
diff --git a/fitpot/read_params.F90 b/fitpot/read_params.F90
index a1a70f8..136b0c5 100644
--- a/fitpot/read_params.F90
+++ b/fitpot/read_params.F90
@@ -299,6 +299,7 @@ subroutine read_params_desc()
   ncnst_type(103) = 1 ! cos(cos(thijk))
   ncnst_type(104) = 1 ! sin(cos(thijk))
   ncnst_type(105) = 2 ! exp(-eta*(cos(thijk)-rs)**2)
+  ncnst_type(106) = 1 ! angular6 (no exp term in angular1)
   ncomb_type(1:100) = 2    ! pair
   ncomb_type(101:200) = 3  ! triplet
 
diff --git a/nappy/io.py b/nappy/io.py
index 77507a2..f17ea5f 100644
--- a/nappy/io.py
+++ b/nappy/io.py
@@ -19,11 +19,11 @@
     scaled_to_cartesian, cartesian_to_scaled
 
 __author__ = "RYO KOBAYASHI"
-__version__ = "240323"
+__version__ = "240903"
 
 READ_FORMATS = ('pmd','POSCAR','CONTCAR','dump','xsf','lammps',
                 'cube','CHGCAR','pdb','extxyz')
-WRITE_FORMATS = ('pmd','POSCAR','dump','xsf','lammps', 'extxyz')
+WRITE_FORMATS = ('pmd','POSCAR','dump','xsf','lammps', 'extxyz', 'pdb')
 
 def write(nsys,fname="pmdini",format=None,**kwargs):
     global myopen, open
@@ -53,8 +53,11 @@ def write(nsys,fname="pmdini",format=None,**kwargs):
         write_cube(nsys,fname,**kwargs)
     elif format in ('pdb','PDB'):
         import ase.io
-        ase.io.write(filename=fname,images=nsys.to_ase_atoms,
+        ase.io.write(filename=fname,images=nsys.to_ase_atoms(),
                      format='proteindatabank')
+    elif format in ('extxyz'):
+        with open(fname,'w') as f:
+            write_extxyz(f,nsys)
     else:
         raise IOError('Cannot write out in the given format: '+format)
 
@@ -90,6 +93,8 @@ def read(fname="pmdini",format=None, specorder=None, index=None):
         nsys = read_lammps_data(fname,specorder=specorder)
     elif format == 'cube':
         nsys = read_cube(fname,specorder=specorder)
+    elif format == 'extxyz':
+        nsys = read_extxyz(fname,specorder=specorder)
     else:
         print('Since the file format is unknown, try to read the file using ASE.')
         try:
@@ -1035,11 +1040,11 @@ def write_extxyz(fileobj, nsys):
     return None
 
 
-def read_extxyz(fileobj, ):
+def read_extxyz(fname, specorder=None):
     """
-    Read a nsys from extxyz format fileobj.
-    Since the extxyz can contain multiple configurations, this method uses
-    fileobj instead of filename.
+    Read an extxyz format using ASE package.
+    NOTE: extxyz file could contain multiple structures.
+    
     The extxyz format is defined in ASE and is like following:
     ---
     8
@@ -1055,33 +1060,21 @@ def read_extxyz(fileobj, ):
     ---
     """
 
-    natm = int(fileobj.readline().split()[0])
-
-    comment = fileobj.readline()
-    # Analyze the comment line
-    
-    
-    hmat = nsys.get_hmat()
-    fileobj.write('Lattice="{0:.3f} {1:.3f} {2:.3f}'.format(*hmat[0,:]))
-    fileobj.write(' {0:.3f} {1:.3f} {2:.3f}'.format(*hmat[1,:]))
-    fileobj.write(' {0:.3f} {1:.3f} {2:.3f}" '.format(*hmat[2,:]))
-    fileobj.write('Properties=species:S:1:pos:R:3:frc:R:3')
-    fileobj.write('\n')
-
-    symbols = nsys.get_symbols()
-    poss = nsys.get_real_positions()
-    frcs = nsys.get_real_forces()
-    for i in range(len(nsys)):
-        si = symbols[i]
-        pi = poss[i]
-        fi = frcs[i]
-        fileobj.write(f'{si:2s}')
-        fileobj.write(f'{pi[0]:16.4f} {pi[1]:16.4f} {pi[2]:16.4f}')
-        fileobj.write(f'{fi[0]:16.3e} {fi[1]:16.3e} {fi[2]:16.3e}')
-        fileobj.write('\n')
+    try:
+        import ase.io
+        atoms = ase.io.read(fname,format='extxyz',index=0)
+        nsys = from_ase(atoms)
+    except Exception as e:
+        print(' Failed to load input file even with ase.')
+        raise
 
-    return None
+    if specorder != None:
+        nsys.specorder != specorder
+        print(' specorder given     = ',specorder)
+        print(' specorder from file = ',nsys.specorder)
+        raise ValueError('Specorder specifically given and obtained from the file do not match!')
 
+    return nsys
     
 def read_CHGCAR(fname='CHGCAR',specorder=None,):
     """
@@ -1280,6 +1273,7 @@ def write_cube(nsys, fname='cube', origin=[0.,0.,0.]):
         f.write(txt)
     return None
 
+
 def get_cube_txt(nsys,origin=[0.,0.,0.]):
     """
     Conver the system info to Gaussian cube format for the purpose of visualization using py3Dmol.
diff --git a/nappy/napsys.py b/nappy/napsys.py
index f70fa1f..ace40f9 100755
--- a/nappy/napsys.py
+++ b/nappy/napsys.py
@@ -420,7 +420,7 @@ def set_real_positions(self,rposs):
         return None
 
     def get_scaled_positions(self):
-        return self.atoms[['x','y','z']].to_numpy()
+        return self.atoms[['x','y','z']].to_numpy(dtype=float)
 
     def set_scaled_positions(self,sposs):
         if len(sposs) != len(self.atoms):
@@ -431,7 +431,7 @@ def set_scaled_positions(self,sposs):
         return None
 
     def get_scaled_velocities(self):
-        return self.atoms[['vx','vy','vz']].to_numpy()
+        return self.atoms[['vx','vy','vz']].to_numpy(dtype=float)
 
     def get_real_velocities(self):
         hmat = self.get_hmat()
@@ -463,7 +463,7 @@ def set_real_velocities(self,vels):
         return None
 
     def get_scaled_forces(self):
-        return self.atoms[['fx','fy','fz']].to_numpy()
+        return self.atoms[['fx','fy','fz']].to_numpy(dtype=float)
 
     def set_scaled_forces(self,sfrcs):
         assert len(sfrcs) == len(self.atoms), 'Array size inconsistent.'
@@ -750,9 +750,9 @@ def get_distance(self,ia,ja):
             raise ValueError('ia > natms, ia,natms = ',ia,natm)
         if ja > natm:
             raise ValueError('ja > natms, ja,natms = ',ja,natm)
-        xi = self.atoms.loc[ia,['x','y','z']].to_numpy()
-        xj = self.atoms.loc[ja,['x','y','z']].to_numpy()
-        xij = xj-xi -np.round(xj-xi)
+        xi = self.atoms.loc[ia,['x','y','z']].to_numpy(dtype=float)
+        xj = self.atoms.loc[ja,['x','y','z']].to_numpy(dtype=float)
+        xij = xj-xi -np.rint(xj-xi)
         hmat = self.get_hmat()
         rij = np.dot(hmat,xij)
         rij2 = rij[0]**2 +rij[1]**2 +rij[2]**2
@@ -769,9 +769,9 @@ def get_angle(self,i,j,k):
             raise ValueError('j > natms, j,natms = ',j,natm)
         if k > natm:
             raise ValueError('k > natms, k,natms = ',k,natm)
-        xi = self.atoms.loc[i,['x','y','z']].to_numpy()
-        xj = self.atoms.loc[j,['x','y','z']].to_numpy()
-        xk = self.atoms.loc[k,['x','y','z']].to_numpy()
+        xi = self.atoms.loc[i,['x','y','z']].to_numpy(dtype=float)
+        xj = self.atoms.loc[j,['x','y','z']].to_numpy(dtype=float)
+        xk = self.atoms.loc[k,['x','y','z']].to_numpy(dtype=float)
         xij = xj-xi -np.round(xj-xi)
         xik = xk-xi -np.round(xk-xi)
         hmat = self.get_hmat()
@@ -1074,8 +1074,8 @@ def repeat(self,n1o,n2o,n3o,n1m=0,n2m=0,n3m=0):
             newauxs[auxname] = []
         inc = 0
         poss = self.get_scaled_positions()
-        vels = self.atoms[['vx','vy','vz']].to_numpy()
-        frcs = self.atoms[['fx','fy','fz']].to_numpy()
+        vels = self.atoms[['vx','vy','vz']].to_numpy(dtype=float)
+        frcs = self.atoms[['fx','fy','fz']].to_numpy(dtype=float)
         for i1 in range(n1m,n1):
             for i2 in range(n2m,n2):
                 for i3 in range(n3m,n3):
@@ -1156,8 +1156,8 @@ def divide(self,*ds):
             newauxs[auxname] = []
         inc = 0
         poss = self.get_scaled_positions()
-        vels = self.atoms[['vx','vy','vz']].to_numpy()
-        frcs = self.atoms[['fx','fy','fz']].to_numpy()
+        vels = self.atoms[['vx','vy','vz']].to_numpy(dtype=float)
+        frcs = self.atoms[['fx','fy','fz']].to_numpy(dtype=float)
         for ia in range(len(self.atoms)):
             pi = poss[ia]
             survive = True
diff --git a/nappy/nn/make_params_NN.py b/nappy/nn/make_params_NN.py
index 87c018d..c56704f 100644
--- a/nappy/nn/make_params_NN.py
+++ b/nappy/nn/make_params_NN.py
@@ -18,7 +18,7 @@
 import json
 
 _descfname = 'in.params.desc'
-_paramfname = 'in.params.NN2'
+_paramfname = 'in.params.'
 _logfname = 'log.make_params_NN'
 
 #...initial range of parameters
@@ -26,7 +26,8 @@
 _pmax =  20.0
 
 _type_avail = ('Gaussian','cosine','polynomial','Morse','angular',
-               'angular1','angular2','angular3','angular4','angular5')
+               'angular1','angular2','angular3','angular4','angular5',
+               'angular6')
 _type_2body = ('Gaussian','cosine','polynomial','Morse',)
 _type_3body = ('angular','angular1','angular2','angular3','angular4','angular5')
 _type2num = {
@@ -40,6 +41,7 @@
     'angular3': 103,
     'angular4': 104,
     'angular5': 105,
+    'angular6': 106,
 }
 _nparam_type = {
     'Gaussian': 2,
@@ -52,6 +54,7 @@
     'angular3': 2,
     'angular4': 2,
     'angular5': 3,
+    'angular6': 3,
 }
 
 def _num2type(num):
@@ -275,7 +278,7 @@ def get_nsf3(triplets):
         nsf3 += len(triplet.sfparams)
     return nsf3
     
-def create_param_files(inputs,nsf2,pairs,nsf3,triplets):
+def create_param_files(inputs,nsf2,pairs,nsf3,triplets,frctype='NN2'):
 
     # print('inputs = ',inputs)
     nl = inputs['nl']
@@ -339,7 +342,7 @@ def create_param_files(inputs,nsf2,pairs,nsf3,triplets):
                 cspi,cspj,cspk = triple.get_csps()
                 for isf,sf in enumerate(triple.sfparams):
                     t = sf[0]
-                    if t in ('angular','angular1'):
+                    if t in ('angular','angular1','angular6'):
                         a = sf[1]
                         ang = -math.cos(a/180*math.pi)
                         f.write(' {0:3d} '.format(_type2num[t]) \
@@ -370,17 +373,27 @@ def create_param_files(inputs,nsf2,pairs,nsf3,triplets):
                                 +' {0:10.4f} {1:10.4f}\n'.format(eta,rs))
         f.close()
 
-    with open(_paramfname,'w') as g:
-        if nl == 1:
-            nc= nhl[0]*nhl[1] +nhl[1]
-            g.write(' {0:4d} {1:6d} {2:4d}\n'.format(nl,nsf,nhl[1]))
-        elif nl == 2:
-            nc= nhl[0]*nhl[1] +nhl[1]*nhl[2] +nhl[2]
-            g.write(' {0:4d} {1:4d} {2:3d} {3:3d}\n'.format(nl,nsf,nhl[1],nhl[2]))
-        for ic in range(nc):
-            g.write(' {0:10.6f}'.format(random.uniform(_pmin,_pmax)))
-            g.write(' {0:10.4f} {1:10.4f}\n'.format(_pmin,_pmax))
-        g.close()
+    if frctype == 'NN2':
+        with open(_paramfname+frctype,'w') as g:
+            if nl == 1:
+                nc= nhl[0]*nhl[1] +nhl[1]
+                g.write(' {0:4d} {1:6d} {2:4d}\n'.format(nl,nsf,nhl[1]))
+            elif nl == 2:
+                nc= nhl[0]*nhl[1] +nhl[1]*nhl[2] +nhl[2]
+                g.write(' {0:4d} {1:4d} {2:3d} {3:3d}\n'.format(nl,nsf,nhl[1],nhl[2]))
+            for ic in range(nc):
+                g.write(' {0:10.6f}'.format(random.uniform(_pmin,_pmax)))
+                g.write(' {0:10.4f} {1:10.4f}\n'.format(_pmin,_pmax))
+            g.close()
+    elif frctype == 'linreg':
+        with open(_paramfname+frctype,'w') as g:
+            nc = nhl[0]
+            g.write(f' {nc:d}\n')
+            for ic in range(nc):
+                g.write(' {0:10.6f}\n'.format(random.uniform(_pmin,_pmax)))
+            g.close()
+    else:
+        raise NotImplementedError(f'No such frctype: {frctype}')
 
     return
 
diff --git a/pmd/descriptor.F90 b/pmd/descriptor.F90
index f2fbbae..9c9f469 100644
--- a/pmd/descriptor.F90
+++ b/pmd/descriptor.F90
@@ -106,6 +106,7 @@ subroutine init_desc()
     ncnst_type(103) = 1 ! cos(cos(thijk))
     ncnst_type(104) = 1 ! sin(cos(thijk))
     ncnst_type(105) = 2 ! exp(-eta*(cos(thijk)-rs)**2)
+    ncnst_type(106) = 1 ! angular6 (no exp term in angular1)
     ncomb_type(1:100) = 2    ! pair
     ncomb_type(101:200) = 3  ! triplet
 
@@ -698,6 +699,42 @@ subroutine calc_desc_default(namax,natm,nb,nnmax,h,tag,ra,lspr,rc &
               igsf(isf,0,ia) = 1
               igsf(isf,jj,ia) = 1
               igsf(isf,kk,ia) = 1
+!.....iyp==106:  (lmd +cos(thijk))^2/(lmd+1)^2 *fcij*fcik
+!.....   for force_fdesc with less dependent on bond length,
+!.....   but rather on bond angles.
+            else if( ityp.eq.106 ) then
+!.....fcij's should be computed after rcs is determined
+              call get_fc_dfc(dij,is,js,rcut,fcij,dfcij)
+              call get_fc_dfc(dik,is,ks,rcut,fcik,dfcik)
+              almbd= desci%prms(1)
+              t2= (abs(almbd)+1d0)**2
+              driki(1:3)= -rik(1:3)/dik
+              drikk(1:3)= -driki(1:3)
+!.....function value
+              t1= (almbd +cs)**2
+!!$              eta3 = 0.5d0 /rcut**2
+!!$              texp = exp(-eta3*(dij2+dik2))
+              tmp = t1/t2 !*texp
+              gsf(isf,ia)= gsf(isf,ia) +tmp*fcij*fcik
+              gijk = gijk +tmp*fcij*fcik
+!.....derivative
+!!$              dgdij= dfcij *fcik *tmp &
+!!$                   +tmp *(-2d0*eta3*dij) *fcij*fcik 
+!!$              dgdik= fcij *dfcik *tmp &
+!!$                   +tmp *(-2d0*eta3*dik) *fcij*fcik 
+              dgdij= dfcij *fcik *tmp
+              dgdik= fcij *dfcik *tmp
+              dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) &
+                   +dgdij*driji(1:3) +dgdik*driki(1:3)
+              dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +dgdij*drijj(1:3)
+              dgsf(1:3,isf,kk,ia)= dgsf(1:3,isf,kk,ia) +dgdik*drikk(1:3)
+              dgcs= 2d0*(almbd+cs)/t2 *fcij*fcik !*texp 
+              dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +dgcs*dcsdi(1:3)
+              dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +dgcs*dcsdj(1:3)
+              dgsf(1:3,isf,kk,ia)= dgsf(1:3,isf,kk,ia) +dgcs*dcsdk(1:3)
+              igsf(isf,0,ia) = 1
+              igsf(isf,jj,ia) = 1
+              igsf(isf,kk,ia) = 1
             endif
 
           enddo ! isf=1,...
@@ -992,6 +1029,42 @@ subroutine desci_bpsf(ia,namax,natm,nnmax,h,tag,ra,lspr,rc,iprint)
             igsfi(isf,0) = 1
             igsfi(isf,jj) = 1
             igsfi(isf,kk) = 1
+!.....iyp==106:  (lmd +cos(thijk))^2/(lmd+1)^2 *fcij*fcik
+!.....   for force_fdesc with less dependent on bond length,
+!.....   but rather on bond angles.
+          else if( ityp.eq.106 ) then
+!.....fcij's should be computed after rcs is determined
+            call get_fc_dfc(dij,is,js,rcut,fcij,dfcij)
+            call get_fc_dfc(dik,is,ks,rcut,fcik,dfcik)
+            almbd= desci%prms(1)
+            t2= (abs(almbd)+1d0)**2
+            driki(1:3)= -rik(1:3)/dik
+            drikk(1:3)= -driki(1:3)
+!.....function value
+            t1= (almbd +cs)**2
+!!$            eta3 = 0.5d0 /rcut**2
+!!$            texp = exp(-eta3*(dij2+dik2))
+            tmp = t1/t2 !*texp
+            gsfi(isf)= gsfi(isf) +tmp*fcij*fcik
+            gijk = gijk +tmp*fcij*fcik
+!.....derivative
+!!$            dgdij= dfcij *fcik *tmp &
+!!$                 +tmp *(-2d0*eta3*dij) *fcij*fcik 
+!!$            dgdik= fcij *dfcik *tmp &
+!!$                 +tmp *(-2d0*eta3*dik) *fcij*fcik 
+            dgdij= dfcij *fcik *tmp
+            dgdik= fcij *dfcik *tmp
+            dgsfi(1:3,isf,0)= dgsfi(1:3,isf,0) &
+                 +dgdij*driji(1:3) +dgdik*driki(1:3)
+            dgsfi(1:3,isf,jj)= dgsfi(1:3,isf,jj) +dgdij*drijj(1:3)
+            dgsfi(1:3,isf,kk)= dgsfi(1:3,isf,kk) +dgdik*drikk(1:3)
+            dgcs= 2d0*(almbd+cs)/t2 *fcij*fcik !*texp 
+            dgsfi(1:3,isf,0)= dgsfi(1:3,isf,0) +dgcs*dcsdi(1:3)
+            dgsfi(1:3,isf,jj)= dgsfi(1:3,isf,jj) +dgcs*dcsdj(1:3)
+            dgsfi(1:3,isf,kk)= dgsfi(1:3,isf,kk) +dgcs*dcsdk(1:3)
+            igsfi(isf,0) = 1
+            igsfi(isf,jj) = 1
+            igsfi(isf,kk) = 1
           endif
 
         enddo ! isf=1,...
diff --git a/pmd/force_UF3.F90 b/pmd/force_UF3.F90
new file mode 100644
index 0000000..7762005
--- /dev/null
+++ b/pmd/force_UF3.F90
@@ -0,0 +1,636 @@
+module UF3
+!-----------------------------------------------------------------------
+!                     Last modified: <2024-10-06 00:43:51 KOBAYASHI Ryo>
+!-----------------------------------------------------------------------
+!  Parallel implementation of Ultra-Fast Force-Field (UF3) for pmd
+!    - 2024.09.02 by R.K., start to implement
+!-----------------------------------------------------------------------
+  use pmdvars,only: nspmax
+  use util,only: csp2isp, num_data
+  use memory,only: accum_mem
+  use vector,only: dot
+  implicit none
+  include "mpif.h"
+  include './const.h'
+  include "./params_unit.h"
+  save
+
+  character(len=128):: paramsdir = '.'
+!.....parameter file name
+  character(128),parameter:: cpfname= 'in.params.uf3'
+  character(128),parameter:: ccfname='in.const.uf3'
+  integer,parameter:: ioprms = 50
+
+  logical:: lprmset_uf3 = .false.
+
+!.....uf2 parameters
+  type prm2
+    character(2):: cb, csi, csj, cknot
+!  cb: NA, 2B or 3B
+!  csi,csj,csk: species name
+!  cknot: nk (non-uniform knot spacing) or uk (uniform knot spacing)
+    integer:: nklead, nktrail
+    integer:: nknot, ncoef
+    real(8):: rc,rc2
+    real(8),allocatable:: knots(:), coefs(:)
+  end type prm2
+
+!.....uf3 parameters
+  type prm3
+    character(2):: cb, csi, csj, csk, cknot
+!  cb: NA, 2B or 3B
+!  csi,csj,csk: species name
+!  cknot: nk (non-uniform knot spacing) or uk (uniform knot spacing)
+    integer:: nklead, nktrail
+    integer:: nknij, nknik, nknjk, ncfij, ncfik, ncfjk
+    real(8):: rcij, rcik, rcjk, rcij2, rcik2, rcjk2
+    real(8),allocatable:: knij(:), knik(:), knjk(:), coefs(:,:,:)
+  end type prm3
+
+  integer:: n2b, n3b
+  type(prm2),allocatable:: prm2s(:)
+  type(prm3),allocatable:: prm3s(:)
+  logical:: has_trios = .false.
+  
+!.....Map of pairs (trios) to parameter set id
+  integer:: interact2(nspmax,nspmax), interact3(nspmax,nspmax,nspmax)
+  
+!.....constants
+  integer:: nelem,nexp,nsp
+
+contains
+  subroutine read_params_uf3(myid,mpi_world,iprint)
+!
+!  Read parameters of uf3 potential from in.params.uf3 file that is given
+!  by uf3/lammps_plugin/scripts/generate_uf3_lammps_pots.py, 
+!-----------------------------------------------------------------------
+!  #UF3 POT UNITS: metal DATE: 2024-09-12 17:12:36 AUTHOR: RK CITATION:
+!  2B W W 0 3 nk
+!  5.5  22
+!  0.001 0.001 0.001 0.001 0.3676 0.7342 ... 5.1334 5.5 5.5 5.5 5.5
+!  18
+!  36.82 36.82 26.69 ... -0.032 0 0 0
+!  #
+!  #UF3 POT UNITS: metal DATE: 2024-09-12 17:12:36 AUTHOR: RK CITATION:
+!  3B W W W 0 3 nk
+!  7.0 3.5 3.5 19 13 13
+!  1.5 1.5 1.5 1.96 2.42 ... 6.54 7 7 7
+!  1.5 1.5 1.5 1.83 2.17 ... 3.17 3.5 3.5 3.5
+!  1.5 1.5 1.5 1.83 2.17 ... 3.17 3.5 3.5 3.5
+!  9 9 15
+!  0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
+!  ...
+!  #
+!-----------------------------------------------------------------------
+!  Each pair or trio is sandwitched with lines beginning with "#".
+!  For trios, in the case of knot information the order is jk,ik,ij,
+!  but for coefficient information the order is ij,ik,jk.
+!  And there is a limitation in the current lammps implementation version
+!  such that r_max_jk = 2*r_max_ij = 2_r_max_ik.
+!
+    implicit none
+    integer,intent(in):: myid,mpi_world,iprint
+    
+    integer:: itmp,ierr,i,j,i2b,i3b
+    integer:: nklead, nktrail
+!  nklead, nktrail: num of leading or trailing knots
+    logical:: lexist
+    character:: fname*128, cmode*4, cb*2, csi*2, csj*2, csk*2, &
+         cknot*2, ctmp*128, cline*128
+!  cmode: none or read
+    
+    if( myid == 0 ) then
+      if( iprint >= ipl_basic ) print '(/,a)',' Read UF3 parameters...'
+      fname = trim(paramsdir)//'/'//trim(cpfname)
+!.....read parameters at the 1st call
+      inquire(file=trim(fname),exist=lexist)
+      if( .not. lexist ) then
+        write(6,'(a)') '   [Error] '//trim(fname)//' does not exist !!!.'
+!!$        call mpi_finalize(ierr)
+        stop
+      endif
+      cmode = 'none'
+      n2b = 0
+      n3b = 0
+      interact2(:,:) = -1
+      interact3(:,:,:) = -1
+      open(ioprms,file=trim(fname),status='old')
+!.....1st, count number of 2B and 3B entries to allocate type objects
+      do while(.true.)
+        read(ioprms,'(a)',end=10) cline
+        if( cline(1:1) == '#' ) then
+          if( cline(1:3) == '#UF' ) then
+            cmode = 'read'
+            cb = 'NA'
+          else
+            cmode = 'none'
+          endif
+          cycle
+        endif
+        if( trim(cmode) == 'none' ) cycle
+        if( cb == 'NA' ) then  ! if cb is not read yet
+          read(cline,*,iostat=ierr) cb
+          if( cb == '2B' ) then
+            n2b = n2b + 1
+          else if( cb == '3B' ) then
+            n3b = n3b + 1
+            has_trios = .true.
+          endif
+! do nothing if cb is already set...
+        endif
+!        if( ierr /= 0 ) cycle
+      end do  ! finished counting 2B & 3B entries
+10    continue
+      if( iprint >= ipl_debug ) then
+        print *,'  n2b, n3b = ',n2b,n3b
+        print *,'  has_trios= ',has_trios
+      endif
+!.....allocate lists of uf2prms and uf3prms
+      if( allocated(prm2s) ) deallocate(prm2s)
+      if( allocated(prm3s) ) deallocate(prm3s)
+      if( n2b /= 0 ) allocate(prm2s(n2b))
+      if( n3b /= 0 ) allocate(prm3s(n3b))
+      rewind(ioprms)
+      i2b = 0
+      i3b = 0
+      do while(.true.)
+        read(ioprms,'(a)',end=20) cline
+        if( cline(1:1) == '#') then
+          if( cline(1:3) == '#UF') then
+            cmode = 'read'
+            cb = 'NA'
+          else
+            cmode = 'none'
+          endif
+          cycle
+        endif
+        if( trim(cmode) == 'none' ) cycle
+        if( cb == 'NA' ) then  ! if cb is not read yet
+          read(cline,*,iostat=ierr) cb
+          if( cb == '2B' ) then
+            backspace(ioprms)
+            i2b = i2b +1
+            call read_2b(prm2s(i2b),i2b)
+            if( iprint >= ipl_debug ) call print_2b(prm2s(i2b))
+          else if( cb == '3B' ) then
+            backspace(ioprms)
+            i3b = i3b +1
+            call read_3b(prm3s(i3b),i3b)
+            if( iprint >= ipl_debug ) call print_3b(prm3s(i3b))
+          endif
+        endif
+!        if( ierr /= 0 ) cycle
+
+      enddo  ! while(.true.)
+20    continue  ! when the file reached the end
+    endif
+
+    call bcast_uf3_params(mpi_world,myid)
+
+    return
+  end subroutine read_params_uf3
+!=======================================================================
+  subroutine read_2b(ps,i2b)
+!
+!  Read 2B part of file via ioprms.
+!  Assuming that the starting line of IOPRMS is the 1st line after
+!  the line starting with '#UF3'.
+!
+    type(prm2),intent(out):: ps
+    integer,intent(in):: i2b
+    integer:: i, isp, jsp
+    
+    read(ioprms,*) ps%cb, ps%csi, ps%csj, ps%nklead, ps%nktrail, ps%cknot
+    if( ps%cb /= '2B' ) stop 'ERROR: CB should be 2B.'
+    read(ioprms,*) ps%rc, ps%nknot
+    ps%rc2 = ps%rc**2
+    if( allocated(ps%knots) ) deallocate(ps%knots)
+    allocate(ps%knots(ps%nknot))
+    read(ioprms,*) (ps%knots(i), i=1,ps%nknot)
+    read(ioprms,*) ps%ncoef
+    if( allocated(ps%coefs) ) deallocate(ps%coefs)
+    allocate(ps%coefs(ps%ncoef))
+    read(ioprms,*) (ps%coefs(i), i=1,ps%ncoef)
+
+    isp = csp2isp(ps%csi)
+    jsp = csp2isp(ps%csj)
+    interact2(isp,jsp) = i2b
+    interact2(jsp,isp) = i2b
+  end subroutine read_2b
+!=======================================================================
+  subroutine read_3b(ps,i3b)
+!
+!  Read 3B part of file via ioprms
+!  Assuming that the starting line of IOPRMS is the 1st line after
+!  the line starting with '#UF3'.
+!
+    type(prm3),intent(out):: ps
+    integer,intent(in):: i3b
+    integer:: i,j,k,isp,jsp,ksp
+    
+    read(ioprms,*) ps%cb, ps%csi, ps%csj, ps%csk, ps%nklead, ps%nktrail, ps%cknot
+    if( ps%cb /= '3B' ) stop 'ERROR: CB should be 3B.'
+    read(ioprms,*) ps%rcjk, ps%rcij, ps%rcik, ps%nknjk, ps%nknij, ps%nknik
+    ps%rcjk2 = ps%rcjk**2
+    ps%rcik2 = ps%rcik**2
+    ps%rcij2 = ps%rcij**2
+    if( allocated(ps%knij) ) deallocate(ps%knij, ps%knik, ps%knjk)
+    allocate(ps%knij(ps%nknij), ps%knik(ps%nknik), ps%knjk(ps%nknjk))
+    read(ioprms,*) (ps%knjk(i), i=1,ps%nknjk)
+    read(ioprms,*) (ps%knij(i), i=1,ps%nknij)
+    read(ioprms,*) (ps%knik(i), i=1,ps%nknik)
+    read(ioprms,*) ps%ncfij, ps%ncfik, ps%ncfjk
+    if( allocated(ps%coefs) ) deallocate(ps%coefs)
+    allocate(ps%coefs(ps%ncfij, ps%ncfik, ps%ncfjk))
+    do i=1, ps%ncfij
+      do j=1, ps%ncfik
+        read(ioprms,*) (ps%coefs(i,j,k), k=1,ps%ncfjk)
+      enddo
+    enddo
+
+    isp = csp2isp(ps%csi)
+    jsp = csp2isp(ps%csj)
+    ksp = csp2isp(ps%csk)
+    interact3(isp,jsp,ksp) = i3b
+    interact3(isp,ksp,jsp) = i3b
+  end subroutine read_3b
+!=======================================================================
+  subroutine bcast_uf3_params(mpi_world,myid)
+!
+!  Broadcast 2B & 3B parameters.
+!
+    integer,intent(in):: mpi_world, myid
+    integer:: i, i2b, i3b, ierr
+    type(prm2):: p2
+    type(prm3):: p3
+
+    call mpi_bcast(n2b, 1, mpi_integer, 0,mpi_world,ierr)
+    call mpi_bcast(n3b, 1, mpi_integer, 0,mpi_world,ierr)
+    if( .not. allocated(prm2s) ) allocate(prm2s(n2b))
+    if( .not. allocated(prm3s) ) allocate(prm3s(n3b))
+    
+    do i2b=1,n2b
+      p2 = prm2s(i2b)
+      call mpi_bcast(p2%cb,2,mpi_character,0,mpi_world,ierr)
+      call mpi_bcast(p2%csi,2,mpi_character,0,mpi_world,ierr)
+      call mpi_bcast(p2%csj,2,mpi_character,0,mpi_world,ierr)
+      call mpi_bcast(p2%cknot,2,mpi_character,0,mpi_world,ierr)
+      call mpi_bcast(p2%nklead,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p2%nktrail,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p2%nknot,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p2%ncoef,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p2%rc,1,mpi_real8,0,mpi_world,ierr)
+      if( .not.allocated(p2%knots) ) allocate(p2%knots(p2%nknot))
+      if( .not.allocated(p2%coefs) ) allocate(p2%coefs(p2%ncoef))
+      call mpi_bcast(p2%knots,p2%nknot,mpi_real8,0,mpi_world,ierr)
+      call mpi_bcast(p2%coefs,p2%ncoef,mpi_real8,0,mpi_world,ierr)
+    enddo
+
+    do i3b=1,n3b
+      p3 = prm3s(i3b)
+      call mpi_bcast(p3%cb,2,mpi_character,0,mpi_world,ierr)
+      call mpi_bcast(p3%csi,2,mpi_character,0,mpi_world,ierr)
+      call mpi_bcast(p3%csj,2,mpi_character,0,mpi_world,ierr)
+      call mpi_bcast(p3%csk,2,mpi_character,0,mpi_world,ierr)
+      call mpi_bcast(p3%cknot,2,mpi_character,0,mpi_world,ierr)
+      call mpi_bcast(p3%nklead,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p3%nktrail,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p3%nknij,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p3%nknik,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p3%nknjk,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p3%ncfij,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p3%ncfik,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p3%ncfjk,1,mpi_integer,0,mpi_world,ierr)
+      call mpi_bcast(p3%rcij,1,mpi_real8,0,mpi_world,ierr)
+      call mpi_bcast(p3%rcik,1,mpi_real8,0,mpi_world,ierr)
+      call mpi_bcast(p3%rcjk,1,mpi_real8,0,mpi_world,ierr)
+      if( .not.allocated(p3%knij) ) allocate(p3%knij(p3%nknij))
+      if( .not.allocated(p3%knik) ) allocate(p3%knik(p3%nknik))
+      if( .not.allocated(p3%knjk) ) allocate(p3%knjk(p3%nknjk))
+      if( .not.allocated(p3%coefs)) allocate(p3%coefs(p3%ncfij, p3%ncfik, p3%ncfjk))
+      call mpi_bcast(p3%knij,p3%nknij,mpi_real8,0,mpi_world,ierr)
+      call mpi_bcast(p3%knik,p3%nknik,mpi_real8,0,mpi_world,ierr)
+      call mpi_bcast(p3%knjk,p3%nknjk,mpi_real8,0,mpi_world,ierr)
+      call mpi_bcast(p3%coefs,p3%ncfjk*p3%ncfik*p3%ncfij,mpi_real8,0,mpi_world,ierr)
+    enddo
+
+    call mpi_bcast(has_trios, 1, mpi_logical, 0,mpi_world,ierr)
+    call mpi_bcast(interact2, nspmax*nspmax, mpi_integer, 0,mpi_world,ierr)
+    call mpi_bcast(interact3, nspmax**3, mpi_integer, 0,mpi_world,ierr)
+    
+  end subroutine bcast_uf3_params
+!=======================================================================
+  subroutine force_uf3(namax,natm,tag,ra,nnmax,aa,strs,h,hi &
+       ,nb,nbmax,lsb,nex,lsrc,myparity,nn,sv,rcin,lspr &
+       ,mpi_world,myid,epi,epot,lstrs,iprint,l1st)
+    use util, only: itotOf
+    implicit none
+    integer,intent(in):: namax,natm,nnmax,iprint
+    integer,intent(in):: nb,nbmax,lsb(0:nbmax,6),lsrc(6),myparity(3) &
+         ,nn(6),mpi_world,myid,lspr(0:nnmax,namax),nex(3)
+    real(8),intent(in):: ra(3,namax),tag(namax) &
+         ,h(3,3),hi(3,3),sv(3,6)
+    real(8),intent(inout):: rcin
+    real(8),intent(out):: aa(3,namax),epi(namax),epot,strs(3,3,namax)
+    logical,intent(in):: l1st 
+    logical:: lstrs
+
+!.....local
+    integer:: ia,ja,ka,jj,kk,l,is,nr2,n,nij3,inc,nik3,njk3,&
+         nik,njk,nij,itot,jtot,ktot,i2b,i3b,js,ks,jsp,ksp,ierr
+    real(8):: epotl2,epotl3,epot2,epot3,tmp,bij,dbij,&
+         bij3(4),dbij3(4),bik3,dbik3,bjk3,dbjk3,c2t,c3t
+    real(8):: xi(3),xj(3),xk(3),xij(3),xik(3),xjk(3),rij(3),rik(3),&
+         rjk(3),dij2,dij,dik2,dik,djk2,djk,drijj(3),drikk(3),&
+         drjkk(3)
+    real(8),save,allocatable:: aal2(:,:),aal3(:,:),strsl(:,:,:)
+    real(8),save:: rcin2
+
+    type(prm2):: p2
+    type(prm3):: p3
+
+    if( l1st ) then
+      if( allocated(aal2) ) deallocate(aal2,strsl)
+      allocate(aal2(3,namax),aal3(3,namax),strsl(3,3,namax))
+      rcin2 = rcin*rcin
+    endif
+
+    if( size(aal2) < 3*namax ) then
+      deallocate(aal2,aal3,strsl)
+      allocate(aal2(3,namax),aal3(3,namax),strsl(3,3,namax))
+    endif
+
+    aal2(:,:) = 0d0
+    aal3(:,:) = 0d0
+    strsl(1:3,1:3,1:namax) = 0d0
+    epotl2 = 0d0
+    epotl3 = 0d0
+    do ia=1,natm
+      is = int(tag(ia))
+      xi(1:3) = ra(1:3,ia)
+      do jj=1,lspr(0,ia)
+        ja = lspr(jj,ia)
+!!$        if( ja <= ia ) cycle
+        js = int(tag(ja))
+!.....Pair terms
+        i2b = interact2(is,js)
+        if( i2b <= 0 ) cycle
+        p2 = prm2s(i2b)
+        xj(1:3) = ra(1:3,ja)
+        xij(1:3) = xj(1:3) -xi(1:3)
+        rij(1:3) = h(1:3,1)*xij(1) +h(1:3,2)*xij(2) +h(1:3,3)*xij(3)
+        dij2 = rij(1)*rij(1) +rij(2)*rij(2) +rij(3)*rij(3)
+        if( dij2 > p2%rc2 ) cycle
+        dij = sqrt(dij2)
+        drijj(1:3) = rij(1:3)/dij
+!.....Look for nr2 from knots data
+        nr2 = knot_index(dij, p2%nknot, p2%knots)
+        do n=nr2-3,nr2
+          if( n < 1 .or. n > p2%nknot-4 ) cycle
+          bij = b_spl(n,3,dij,p2%knots,p2%nknot)
+          c2t = p2%coefs(n)
+          tmp = c2t *bij
+          epi(ia) = epi(ia) +tmp
+          epotl2 = epotl2 + tmp
+!!$          epi(ja) = epi(ja) +tmp
+!!$          epotl2 = epotl2 + tmp + tmp
+!.....Forces
+          dbij = db_spl(n,3,dij,p2%knots,p2%nknot)
+          aal2(1:3,ia) = aal2(1:3,ia) +drijj(1:3)*dbij*c2t
+          aal2(1:3,ja) = aal2(1:3,ja) -drijj(1:3)*dbij*c2t
+        enddo
+      enddo
+
+!.....Trio part is separated from pair part,
+!.....which may be slower because of double computation of dij,
+!.....but this code is a bit simpler.
+      if( .not.has_trios ) cycle
+      do jsp=1,nspmax
+        do ksp=jsp,nspmax
+          i3b = interact3(is,jsp,ksp)
+          if( i3b <= 0 ) cycle
+          p3 = prm3s(i3b)
+          do jj=1,lspr(0,ia)
+            ja = lspr(jj,ia)
+            js = int(tag(ja))
+            if( js /= jsp ) cycle
+            xj(1:3) = ra(1:3,ja)
+            xij(1:3) = xj(1:3) -xi(1:3)
+            rij(1:3) = h(1:3,1)*xij(1) +h(1:3,2)*xij(2) +h(1:3,3)*xij(3)
+            dij2 = rij(1)*rij(1) +rij(2)*rij(2) +rij(3)*rij(3)
+            if( dij2 > p3%rcij2 ) cycle
+            dij = sqrt(dij2)
+            drijj(1:3) = rij(1:3)/dij
+            nij3 = knot_index(dij, p3%nknij, p3%knij)
+            inc = 0
+            bij3(:) = 0d0
+            dbij3(:) = 0d0
+            do n=nij3-3,nij3
+              inc = inc +1
+              if( n < 1 .or. n > p3%nknij-4 ) cycle
+              bij3(inc) = b_spl(n,3,dij, p3%knij, p3%nknij)
+              dbij3(inc) = db_spl(n,3,dij, p3%knij, p3%nknij)
+            enddo
+
+            do kk=1,lspr(0,ia)
+              ka = lspr(kk,ia)
+              if( jsp == ksp .and. ka <= ja ) cycle
+              ks = int(tag(ka))
+              if( ks /= ksp ) cycle
+              xk(1:3) = ra(1:3,ka)
+              xik(1:3) = xk(1:3) -xi(1:3)
+              xjk(1:3) = xk(1:3) -xj(1:3)
+              rik(1:3) = h(1:3,1)*xik(1) +h(1:3,2)*xik(2) +h(1:3,3)*xik(3)
+              dik2 = rik(1)*rik(1) +rik(2)*rik(2) +rik(3)*rik(3)
+              if( dik2 > p3%rcik2 ) cycle
+              dik = sqrt(dik2)
+              rjk(1:3) = h(1:3,1)*xjk(1) +h(1:3,2)*xjk(2) +h(1:3,3)*xjk(3)
+              djk2 = rjk(1)*rjk(1) +rjk(2)*rjk(2) +rjk(3)*rjk(3)
+              if( djk2 > p3%rcjk2 ) cycle
+              djk = sqrt(djk2)
+              drikk(1:3) = rik(1:3)/dik
+              drjkk(1:3) = rjk(1:3)/djk
+              nik3 = knot_index(dik, p3%nknik, p3%knik)
+              njk3 = knot_index(djk, p3%nknjk, p3%knjk)
+!.....B-spline part
+              do nik=nik3-3,nik3
+                if( nik < 1 .or. nik > p3%nknik-4 ) cycle
+                bik3 = b_spl(nik,3,dik, p3%knik, p3%nknik)
+                dbik3 = db_spl(nik,3,dik, p3%knik, p3%nknik)
+                do njk=njk3-3,njk3
+                  if( njk < 1 .or. njk > p3%nknjk-4 ) cycle
+                  bjk3 = b_spl(njk,3,djk, p3%knjk, p3%nknjk)
+                  dbjk3 = db_spl(njk,3,djk, p3%knjk, p3%nknjk)
+                  l = 0
+                  do nij=nij3-3,nij3
+                    l = l +1
+                    if( nij < 1 .or. nij > p3%nknij-4 ) cycle
+!.....Energy
+                    c3t = p3%coefs(nij,nik,njk)
+                    tmp = c3t*bij3(l)*bik3*bjk3
+                    epi(ia) = epi(ia) +tmp
+                    epotl3 = epotl3 +tmp
+!.....Force
+                    aal3(1:3,ia)= aal3(1:3,ia) &
+                         +c3t*(dbij3(l)*bik3*bjk3*drijj(1:3) &
+                         +bij3(l)*dbik3*bjk3*drikk(1:3))
+                    aal3(1:3,ja)= aal3(1:3,ja) &
+                         +c3t*(dbij3(l)*bik3*bjk3*(-drijj(1:3)) &
+                         +bij3(l)*bik3*dbjk3*drjkk(1:3))
+                    aal3(1:3,ka)= aal3(1:3,ka) &
+                         +c3t*(bij3(l)*dbik3*bjk3*(-drikk(1:3)) &
+                         +bij3(l)*bik3*dbjk3*(-drjkk(1:3)))
+                  enddo ! nij
+                enddo ! njk
+              enddo ! nik
+
+            enddo  ! kk
+          enddo  ! jj
+        enddo  ! ksp
+      enddo  ! jsp
+
+    enddo ! ia
+
+    call copy_dba_bk(namax,natm,nbmax,nb,lsb,nex,lsrc,myparity &
+         ,nn,mpi_world,aal2,3)
+    if( has_trios ) call copy_dba_bk(namax,natm,nbmax,nb,lsb,nex,lsrc,myparity &
+         ,nn,mpi_world,aal3,3)
+    aa(1:3,1:natm) = aa(1:3,1:natm) +aal2(1:3,1:natm) +aal3(1:3,1:natm)
+    call copy_dba_bk(namax,natm,nbmax,nb,lsb,nex,lsrc,myparity &
+         ,nn,mpi_world,strsl,9)
+    strs(1:3,1:3,1:natm) = strs(1:3,1:3,1:natm) +strsl(1:3,1:3,1:natm)*0.5d0
+
+!-----gather epot
+    epot2 = 0d0
+    epot3 = 0d0
+    call mpi_allreduce(epotl2,epot2,1,mpi_real8,mpi_sum,mpi_world,ierr)
+    if( has_trios ) call mpi_allreduce(epotl3,epot3,1,mpi_real8,mpi_sum,mpi_world,ierr)
+    epot= epot +epot2 +epot3
+    if( myid == 0 .and. iprint > 2 ) &
+         print '(a,2es12.4)',' force_uf3 epot2,epot3 = ',epot2,epot3
+
+    return
+  end subroutine force_uf3
+!=======================================================================
+  recursive function b_spl(n,d,r,ts,nmax) result(val)
+!
+!  TODO: check the efficiency of recursive func
+!  
+!  Recursive implementation of B-spline function with N and D as indices
+!  and R as an argument.
+!  TS --- a list of {t_n}
+!
+    integer,intent(in):: n,d,nmax
+    real(8),intent(in):: r,ts(nmax)
+    real(8):: val
+    real(8):: denom1,denom2
+    
+    if( d == 0 ) then
+      if( r >= ts(n) .and. r < ts(n+1) ) then
+        val = 1d0
+        return
+      else
+        val = 0d0
+        return
+      endif
+    else
+      val = 0d0
+      denom1 = ts(n+d)-ts(n)
+      if( abs(denom1) > 1d-8 ) then
+        val = val +(r-ts(n))/denom1 *b_spl(n,d-1,r,ts,nmax)
+      endif
+      denom2 = ts(n+d+1)-ts(n+1)
+      if( abs(denom2) > 1d-8 ) then
+        val = val +(ts(n+d+1)-r)/denom2 *b_spl(n+1,d-1,r,ts,nmax)
+      endif
+    endif
+    return
+  end function b_spl
+!=======================================================================
+  function db_spl(n,d,r,ts,nmax)
+    integer,intent(in):: n,d,nmax
+    real(8),intent(in):: r,ts(nmax)
+    real(8):: db_spl
+    real(8):: denom1, denom2
+
+    db_spl = 0d0
+    denom1 = ts(n+d)-ts(n)
+    if( denom1 > 1d-8 ) then
+      db_spl = db_spl +d*b_spl(n,d-1,r,ts,nmax) /denom1
+    endif
+    denom2 = ts(n+d+1)-ts(n+1)
+    if( denom2 > 1d-8 ) then
+      db_spl = db_spl -d*b_spl(n+1,d-1,r,ts,nmax)/denom2
+    endif
+    return
+  end function db_spl
+!=======================================================================
+  function knot_index(r, nknot, knots) result(n)
+    real(8),intent(in):: r
+    integer,intent(in):: nknot
+    real(8),intent(in):: knots(nknot)
+    integer:: n, i
+
+    n = 0
+    do i=1,nknot
+      if( r > knots(i) ) then
+        n = i
+      else
+        return
+      endif
+    enddo
+    return
+  end function knot_index
+!=======================================================================
+  subroutine set_paramsdir_uf3(dname)
+!
+!  Accessor routine to set paramsdir.
+!
+    implicit none
+    character(len=*),intent(in):: dname
+
+    paramsdir = trim(dname)
+    return
+  end subroutine set_paramsdir_uf3
+!=======================================================================
+  subroutine print_2b(ps)
+    type(prm2),intent(in):: ps
+    integer:: i
+    character:: c*2
+    print '(/,a)','   UF3 parameters of 2B for '&
+         //trim(ps%csi)//'-'//trim(ps%csj)
+    print '(5(a,1x))','     cb,csi,csj,cknot = ',ps%cb,ps%csi,ps%csj,ps%cknot
+    print '(a,2i3)','     nklead,nktrail = ',ps%nklead,ps%nktrail
+    print '(a,2i3)','     nknot,ncoef = ',ps%nknot, ps%ncoef
+    print '(a,f6.3)','     rc = ',ps%rc
+    write(c,'(i2)') size(ps%knots)
+    print '(a,'//c//'(1x,f7.3))','     knots =',(ps%knots(i),i=1,size(ps%knots))
+    write(c,'(i2)') size(ps%coefs)
+    print '(a,'//trim(c)//'(1x,f7.3))','     coefs =',(ps%coefs(i),i=1,size(ps%coefs))
+  end subroutine print_2b
+!=======================================================================
+  subroutine print_3b(ps)
+    type(prm3),intent(in):: ps
+    integer:: i
+    character:: c*2
+    print '(/,a)','   UF3 parameters of 3B for '&
+         //trim(ps%csi)//'-'//trim(ps%csj)//'-'//trim(ps%csk)
+    print '(6(a,1x))','     cb,csi,csj,csk,cknot = ',ps%cb,ps%csi,ps%csj,ps%csk,ps%cknot
+    print '(a,2i3)','     nklead,nktrail = ',ps%nklead,ps%nktrail
+    print '(a,6i3)','     nknij,nknik,nknjk,ncfij,ncfik,ncfjk = ', &
+         ps%nknij, ps%nknik, ps%nknjk, ps%ncfij, ps%ncfik, ps%ncfjk
+    print '(a,3f6.3)','     rcij,rcik,rcjk = ',ps%rcij,ps%rcik,ps%rcjk
+    write(c,'(i2)') ps%nknij
+    print '(a,'//c//'(1x,f7.3))','     knij =',(ps%knij(i),i=1,ps%nknij)
+    write(c,'(i2)') ps%nknik
+    print '(a,'//c//'(1x,f7.3))','     knik =',(ps%knik(i),i=1,ps%nknik)
+    write(c,'(i2)') ps%nknjk
+    print '(a,'//c//'(1x,f7.3))','     knjk =',(ps%knjk(i),i=1,ps%nknjk)
+!!$    write(c,'(i2)') size(ps%coefs)
+!!$    print '(a,'//trim(c)//'es11.2)','     coefs =',(ps%coefs(i),i=1,size(ps%coefs))
+  end subroutine print_3b
+end module UF3
+!-----------------------------------------------------------------------
+!     Local Variables:
+!     compile-command: "make pmd"
+!     End:
diff --git a/pmd/force_common.F90 b/pmd/force_common.F90
index 790794b..06ac091 100644
--- a/pmd/force_common.F90
+++ b/pmd/force_common.F90
@@ -46,6 +46,7 @@ subroutine get_force(l1st,epot,stnsr)
   use RFMEAM,only: force_RFMEAM
   use Pellenq,only: force_Pellenq
   use repel,only: force_repel
+  use UF3,only: force_uf3
   use fdesc, only: force_fdesc
   use time, only: accum_time
   implicit none
@@ -441,6 +442,14 @@ subroutine get_force(l1st,epot,stnsr)
     call accum_time('force_Coulomb',mpi_wtime() -tmp)
   endif
 
+  if( use_force('UF3').or.use_force('uf3') ) then
+    tmp = mpi_wtime()
+    call force_uf3(namax,natm,tag,ra,nnmax,aa,strs &
+       ,h,hi,nb,nbmax,lsb,nex,lsrc,myparity,nn,sv,rc,lspr &
+       ,mpi_md_world,myid_md,epi,epot,lstrs,iprint,l1st)
+    call accum_time('force_UF3',mpi_wtime() -tmp)
+  endif
+
   if( use_force('fdesc') ) then
     call force_fdesc(namax,natm,nnmax,lspr,rc,h,hi,tag,ra, &
          aa,epot,aux(iaux_edesc,:),strs, &
@@ -485,6 +494,7 @@ subroutine init_force(linit)
   use angular,only: read_params_angular, lprmset_angular
   use RFMEAM, only: read_params_RFMEAM, lprmset_RFMEAM
   use Pellenq,only: read_params_Pellenq, lprmset_Pellenq
+  use UF3,only: read_params_uf3, lprmset_uf3
   use repel,only: read_params_repel, lprmset_repel
   use fdesc, only: init_fdesc
   implicit none
@@ -665,6 +675,13 @@ subroutine init_force(linit)
     endif
   endif
 
+!.....UF3
+  if( use_force('UF3') .or. use_force('uf3') ) then
+    if( .not.lprmset_uf3 ) then
+      call read_params_uf3(myid_md,mpi_md_world,iprint)
+    endif
+  endif
+
   if( use_force('fdesc') ) then
     call init_fdesc(myid_md,mpi_md_world,iprint)
   endif
diff --git a/pmd/force_fdesc.F90 b/pmd/force_fdesc.F90
index d6cbf85..5ee4ceb 100644
--- a/pmd/force_fdesc.F90
+++ b/pmd/force_fdesc.F90
@@ -1,6 +1,6 @@
 module fdesc
 !-----------------------------------------------------------------------
-!                     Last-modified: <2024-07-28 11:51:17 KOBAYASHI Ryo>
+!                     Last-modified: <2024-08-23 16:13:09 KOBAYASHI Ryo>
 !-----------------------------------------------------------------------
 !  Potential in descriptor space.
 !  Originally for the purpose of restricting structure, at 2021-05-17, by R.K.
@@ -18,14 +18,17 @@ module fdesc
 
   logical:: lfdesc = .false.  ! Flag to use fdesc.
   logical:: initialized = .false.
-  logical:: pc1_as_edesc = .false.  ! pc1 value at edesci
+!!$  logical:: pc1_as_edesc = .false.  ! pc1 value at edesci
+! edesc_type: 0) energy, 1) dsq
+  integer:: edesc_type  = 0 
 
   integer:: ndim_desc = -1
   integer:: giddesc = -1
   real(8):: scnst = 1.0d0
   real(8):: gcoef = 0.1d0
   real(8):: gsgm  = 1.0d0
-  real(8),allocatable:: desctgt(:),descpca(:), descfrc(:,:),descstrs(:,:,:)
+  real(8),allocatable:: desctgt(:), descov(:,:), descacc(:,:), &
+       descpca(:), descfrc(:,:),descstrs(:,:,:)
   real(8):: edesc
 !  logical:: ldspc(nspmax)
 
@@ -56,6 +59,7 @@ subroutine init_fdesc(myid,mpi_world,iprint)
 !        print '(a,es11.3)','   Spring constant      = ',scnst
         print '(a,f7.3)','   Gaussian coeff.      = ',gcoef
         print '(a,f7.3)','   Gaussian sigma       = ',gsgm
+        print '(a,i3)',  '   edesc type           = ',edesc_type
       endif
     endif
 
@@ -81,18 +85,18 @@ subroutine read_fdesc_params(myid,mpi_world,iprint)
 !  ! gauss_coef         1.0   (common for all the descriptors)
 !  ! gauss_sigma        1.0   (common for all the descriptors)
 !  desc_dimension     55
-!  !  isf  desctgt  descpca
-!     1    0.1245   0.345  
-!     2    0.1452   0.456  
-!     3    0.1542   0.567  
+!  !  desctgt: (desctgt(i),i=1,ndim)
+!     0.123  0.234  2.345  3.212  0.432 ...
+!  !  descov: ((descov(i,j),j=1,ndim),i=1,ndim)  ! nsf entries per line, nsf lines
+!     0.123  0.234  0.345  ...
 !     ....
-!     55   0.1245   1.234  
+!     0.987  0.876  0.765 ...
 !-----------------------------------------------------------------------
 !
     use util, only: num_data
     integer,intent(in):: myid,mpi_world,iprint
     
-    integer:: nentry,isp,ierr,isf,idesc
+    integer:: nentry,isp,ierr,isf,jsf,idesc
     real(8):: tmp
     character:: cline*128, c1st*128, csp*3, cmode*128
     real(8),allocatable:: desctmp(:)
@@ -121,26 +125,30 @@ subroutine read_fdesc_params(myid,mpi_world,iprint)
           else if( trim(c1st).eq.'gauss_sigma' ) then
             backspace(ionum)
             read(ionum,*) c1st, gsgm
-          else if( trim(c1st).eq.'pc1_as_edesc' ) then
+          else if( trim(c1st).eq.'edesc_type' ) then
             backspace(ionum)
-            read(ionum,*) c1st, pc1_as_edesc
+            read(ionum,*) c1st, edesc_type
           else if( trim(c1st).eq.'desc_dimension' ) then
             backspace(ionum)
             read(ionum,*) c1st, ndim_desc
-            if( allocated(desctgt) ) deallocate(desctgt,descpca)
+            if( allocated(desctgt) ) deallocate(desctgt,descov,descacc)
             allocate(desctmp(ndim_desc),desctgt(ndim_desc), &
-                 descpca(ndim_desc))
+                 descov(ndim_desc,ndim_desc), descacc(ndim_desc,ndim_desc))
             cmode = 'read_descs'
           else
             print *,'There is no such fdesc keyword: ', trim(c1st)
           endif
         else if( trim(cmode).eq.'read_descs' ) then
           backspace(ionum)
-          read(ionum,*) idesc, desctgt(idesc), descpca(idesc)
-          if( idesc.gt.ndim_desc ) then
-            print *,'ERROR: desc-ID exceeds the descriptor dimension.'
-            stop 1
-          endif
+!.....Read target descriptor values
+          read(ionum,*) (desctgt(isf),isf=1,ndim_desc)
+!.....Read covariance matrix of the descriptor
+          read(ionum,*) ((descov(isf,jsf),jsf=1,ndim_desc),isf=1,ndim_desc)
+!!$          read(ionum,*) idesc, desctgt(idesc), descpca(idesc)
+!!$          if( idesc.gt.ndim_desc ) then
+!!$            print *,'ERROR: desc-ID exceeds the descriptor dimension.'
+!!$            stop 1
+!!$          endif
         endif
       enddo
 !!$      read(ionum,*) tmp
@@ -152,14 +160,18 @@ subroutine read_fdesc_params(myid,mpi_world,iprint)
     call mpi_bcast(ndim_desc,1,mpi_integer,0,mpi_world,ierr)
     call mpi_bcast(scnst,1,mpi_real8,0,mpi_world,ierr)
     call mpi_bcast(giddesc,1,mpi_integer,0,mpi_world,ierr)
-    call mpi_bcast(pc1_as_edesc,1,mpi_logical,0,mpi_world,ierr)
+    call mpi_bcast(edesc_type,1,mpi_integer,0,mpi_world,ierr)
 !    call mpi_bcast(ldspc,nspmax,mpi_logical,0,mpi_world,ierr)
     if( myid.ne.0) then
-      if( allocated(desctgt) ) deallocate(desctgt, descpca)
-      allocate(desctgt(ndim_desc), descpca(ndim_desc))
+      if( allocated(desctgt) ) deallocate(desctgt, descov, descacc)
+      allocate(desctgt(ndim_desc), descov(ndim_desc,ndim_desc), &
+           descacc(ndim_desc,ndim_desc))
     endif
     call mpi_bcast(desctgt,ndim_desc,mpi_real8,0,mpi_world,ierr)
-    call mpi_bcast(descpca,ndim_desc,mpi_real8,0,mpi_world,ierr)
+!!$    call mpi_bcast(descpca,ndim_desc,mpi_real8,0,mpi_world,ierr)
+    call mpi_bcast(descov,ndim_desc**2,mpi_real8,0,mpi_world,ierr)
+!.....Calculate descacc by inverting descov
+    call ludc_inv(ndim_desc, descov, descacc)
 
     if( allocated(desctmp) ) deallocate(desctmp)
     return
@@ -181,9 +193,10 @@ subroutine force_fdesc(namax,natm,nnmax,lspr,rcin,h,hi,tag,ra, &
     real(8),intent(inout):: aa(3,namax),edesci(namax),epot,strs(3,3,namax)
     logical,intent(in):: l1st
 
-    integer:: ia,igv,isf,isp,jsp,jj,ja,i,k,ixyz,jxyz,ierr
+    integer:: ia,igv,isf,jsf,isp,jsp,jj,ja,i,k,ixyz,jxyz,ierr
     real(8):: tmp,at(3),ave_aa,ave_dsp,dsq,dpca1,pca1,aexp, &
-         xi(3),xj(3),xij(3),rij(3),dij,sij
+         xi(3),xj(3),xij(3),rij(3),dij,sij,esp
+    real(8):: dxmah(nsf)
 
     if( .not.allocated(descfrc) ) then
       allocate(descfrc(3,namax),descstrs(3,3,namax))
@@ -204,20 +217,28 @@ subroutine force_fdesc(namax,natm,nnmax,lspr,rcin,h,hi,tag,ra, &
       igv = igvarOf(tag(ia),giddesc)
       if( igv.eq.0 ) cycle  ! fdesc works only on atoms of igv > 0 
       isp = int(tag(ia))
-!      if( .not. ldspc(isp) ) cycle  ! only specified species pass here
+!!$      if( .not. ldspc(isp) ) cycle  ! only specified species pass here
       call calc_desci(ia,namax,natm,nnmax,h,tag,ra,lspr,rcin,iprint)
       xi(1:3) = ra(1:3,ia)
-      dpca1 = 0d0
-      pca1 = 0d0
+!!$      dpca1 = 0d0
+!!$      pca1 = 0d0
+      dsq = 0d0
+      dxmah(:) = 0d0
       do isf=1,nsf
-        dpca1 = dpca1 + (gsfi(isf)-desctgt(isf)) * descpca(isf)
-        pca1 = pca1 + gsfi(isf) * descpca(isf)
+!!$        dpca1 = dpca1 + (gsfi(isf)-desctgt(isf)) * descpca(isf)
+!!$        pca1 = pca1 + gsfi(isf) * descpca(isf)
+        do jsf=1,nsf
+          dxmah(isf) = dxmah(isf) +descacc(jsf,isf) *(gsfi(jsf)-desctgt(jsf))
+          dsq = dsq +(gsfi(isf)-desctgt(isf))*descacc(jsf,isf)*(gsfi(jsf)-desctgt(jsf))
+        enddo
       enddo
-      dsq = dpca1**2
-      aexp = -gcoef * exp(-dsq /2 /gsgm**2)
+!!$      dsq = dpca1**2
+!!$      aexp = -gcoef * exp(-dsq /2 /gsgm**2)
+      esp = 0.5d0 *scnst *dsq
       do isf=1,nsf
 !!$        tmp = scnst*dpca1 * descpca(isf)
-        tmp = -aexp *dpca1 /gsgm**2 *descpca(isf)
+!!$        tmp = -aexp *dpca1 /gsgm**2 *descpca(isf)
+        tmp = scnst *dxmah(isf)
 !.....Compute spring forces
         do jj=1,lspr(0,ia)
           ja = lspr(jj,ia)
@@ -238,13 +259,14 @@ subroutine force_fdesc(namax,natm,nnmax,lspr,rcin,h,hi,tag,ra, &
         enddo ! jj
         descfrc(1:3,ia) = descfrc(1:3,ia) -dgsfi(1:3,isf,0)*tmp
       enddo ! isf
-      if( pc1_as_edesc ) then
-        edesci(ia) = edesci(ia) + dpca1
+      if( edesc_type.eq.1 ) then  ! sqrt(dsq) value
+        edesci(ia) = edesci(ia) + sqrt(dsq)
       else
-        edesci(ia) = edesci(ia) + aexp
+        edesci(ia) = edesci(ia) + esp
       endif
 !!$      edesc = edesc +0.5d0 *scnst *dsq
-      edesc = edesc +aexp
+!!$      edesc = edesc +aexp
+      edesc = edesc +esp
     enddo ! ia
 
 !.....Send back forces on immigrants to the neighboring nodes
@@ -255,12 +277,12 @@ subroutine force_fdesc(namax,natm,nnmax,lspr,rcin,h,hi,tag,ra, &
 !!$      at(1:3)= descfrc(1:3,i)
 !!$      descfrc(1:3,i)= hi(1:3,1)*at(1) +hi(1:3,2)*at(2) +hi(1:3,3)*at(3)
 !!$    enddo
-    if( iprint.ge.ipl_debug .and. myid.eq.0 ) then
-      print '(a,es12.3)','edesc = ',edesc
-      do i=1,10
-        print '(i6, 6(2x, es12.3))', i, descfrc(1:3,i), aa(1:3,i)
-      enddo
-    endif
+!!$    if( iprint.ge.ipl_debug .and. myid.eq.0 ) then
+!!$      print '(a,es12.3)','edesc = ',edesc
+!!$      do i=1,10
+!!$        print '(i6, 6(2x, es12.3))', i, descfrc(1:3,i), aa(1:3,i)
+!!$      enddo
+!!$    endif
 !.....Add desc forces to the original forces
     aa(1:3,1:natm) = aa(1:3,1:natm) +descfrc(1:3,1:natm)
 !.....Gather epot
diff --git a/pmd/makefile.in b/pmd/makefile.in
index 40e5f25..27b229f 100644
--- a/pmd/makefile.in
+++ b/pmd/makefile.in
@@ -47,7 +47,7 @@ forces= force_Ramas_FeH.o force_RK_FeH.o force_RK_WHe.o \
 	force_Morse.o force_Coulomb.o force_EAM.o force_Buckingham.o \
 	force_Bonny_WRe.o force_ZBL.o force_cspline.o force_Tersoff.o \
 	force_Abell.o force_BMH.o force_dipole.o force_fpc.o force_angular.o \
-	force_RFMEAM.o force_Pellenq.o force_repel.o force_fdesc.o
+	force_RFMEAM.o force_Pellenq.o force_repel.o force_fdesc.o force_UF3.o
 params=params_Ramas_FeH.h params_RK_FeH.h params_RK_WHe.h params_Ito3_WHe.h \
 	params_EDIP_Si.h \
 	params_Brenner.h params_Lu_WHe.h params_Branicio_AlN.h \
@@ -151,6 +151,7 @@ force_fpc.o: ${mods}
 force_Morse.o: ${mods} force_Coulomb.o
 force_angular.o: ${mods}
 force_fdesc: mod_pmdvars.o mod_util.o descriptor.o
+force_UF3.o: ${mods}
 mod_constraints.o: mod_vector.o mod_util.o
 mod_rdcfrc.o: mod_structure.o mod_util.o
 mod_ttm.o: mod_pmdvars.o mod_util.o mod_memory.o mod_random.o
diff --git a/pmd/mod_force.F90 b/pmd/mod_force.F90
index 34df9a4..e09e16a 100644
--- a/pmd/mod_force.F90
+++ b/pmd/mod_force.F90
@@ -1,6 +1,6 @@
 module force
 !-----------------------------------------------------------------------
-!                     Last-modified: <2024-07-28 11:47:09 KOBAYASHI Ryo>
+!                     Last-modified: <2024-09-28 11:27:25 KOBAYASHI Ryo>
 !-----------------------------------------------------------------------
   use pmdvars,only: nspmax
   implicit none
@@ -13,7 +13,7 @@ module force
        ol_type, ol_force, ol_ranges, ol_alphas, ol_dalphas, ol_pair
   
 !.....Force index list
-  integer,parameter:: N_FORCES = 35
+  integer,parameter:: N_FORCES = 36
   character(len=128):: force_index_list(N_FORCES)
 
   integer:: num_forces = -1
@@ -71,6 +71,7 @@ subroutine init_mod_force()
     force_index_list(33) =  "angular"
     force_index_list(34) =  "DNN"
     force_index_list(35) =  "fdesc"
+    force_index_list(36) =  "fdesc"
 
     luse_force(:) = .false.
 
@@ -121,7 +122,7 @@ subroutine write_forces(myid)
 
     if( myid.eq.0 ) then
       print *,''
-      write(6,'(a)',advance='no') ' Use the following force-fields:'
+      write(6,'(a)',advance='no') ' Force-fields:'
       do i=1,num_forces
         write(6,'(2x,a)',advance='no') trim(force_list(i))
       enddo
diff --git a/pmd/mod_group.F90 b/pmd/mod_group.F90
index d3ef595..e202952 100644
--- a/pmd/mod_group.F90
+++ b/pmd/mod_group.F90
@@ -1,6 +1,6 @@
 module group
 !-----------------------------------------------------------------------
-!                     Last modified: <2024-07-25 11:14:34 KOBAYASHI Ryo>
+!                     Last modified: <2024-07-28 14:30:15 KOBAYASHI Ryo>
 !-----------------------------------------------------------------------
 !  Module for grouping atoms.
 !-----------------------------------------------------------------------
@@ -102,7 +102,9 @@ subroutine register_group(cline)
       do i=1,npair
         ispc(igrp,itmp(i,1)) = itmp(i,2)
       enddo
-      gtiming(igrp) = -1  ! only at the 1st time
+!.....If gtiming==0, set -1 (only at 1st time),
+!.....otherwise keep the value as is.
+      if( gtiming(igrp).eq.0 ) gtiming(igrp) = -1
       deallocate(itmp)
     else if( trim(cgname).eq.'sphere' ) then
 ! e.g.) group  1  sphere  20.0  0.5  0.5  0.5  1  2
@@ -121,6 +123,7 @@ subroutine register_group(cline)
       sph_org(igrp,1:3) = (/ orgx, orgy, orgz /)
       sph_rad2(igrp) = rad**2
       isph(igrp,1:2) = (/ kin, kout /)
+!.....In any case, set gtiming to 1.
       gtiming(igrp) = 1
     else
       print *,'WARNING: no such group name: '//trim(cgname)
diff --git a/pmd/mod_version.F90 b/pmd/mod_version.F90
index 2125aec..2ecc304 100644
--- a/pmd/mod_version.F90
+++ b/pmd/mod_version.F90
@@ -1,6 +1,6 @@
 module version
 !-----------------------------------------------------------------------
-!                     Last-modified: <2023-04-05 13:42:13 KOBAYASHI Ryo>
+!                     Last-modified: <2024-10-06 00:20:56 KOBAYASHI Ryo>
 !-----------------------------------------------------------------------
 ! A module for version/revision.
 !-----------------------------------------------------------------------
@@ -10,7 +10,7 @@ module version
 
   public:: write_revision, write_authors
   
-  character(len=128),parameter:: crevision = 'rev230405'
+  character(len=128),parameter:: crevision = 'rev241006'
 
   character(len=128),parameter:: cauthors(1) = &
        (/ 'Ryo KOBAYASHI <kobayashi.ryo@nitech.ac.jp>' /)
diff --git a/pmd/pmd_core.F90 b/pmd/pmd_core.F90
index 480fb9e..5896123 100644
--- a/pmd/pmd_core.F90
+++ b/pmd/pmd_core.F90
@@ -1,5 +1,5 @@
 !-----------------------------------------------------------------------
-!                     Last-modified: <2024-07-28 11:48:15 KOBAYASHI Ryo>
+!                     Last-modified: <2024-10-05 14:42:10 KOBAYASHI Ryo>
 !-----------------------------------------------------------------------
 ! Core subroutines/functions needed for pmd.
 !-----------------------------------------------------------------------
@@ -572,6 +572,7 @@ subroutine pmd_core(hunit,hmat,ntot0,tagtot,rtot,vtot,atot,stot &
       is = int(tag(i))
       va(1:3,i)=va(1:3,i) +aa(1:3,i)*fa2v(is)*dt
     enddo
+    
     if( chgopt_method(1:4).eq.'xlag' ) call update_vauxq(aux(iaux_vq,:))
     if( index(cpctl,'lange').ne.0 ) call cvel_update_langevin(stnsr,h,mpi_md_world,1)