implement UF3 force field and some bug fixes.

fitpot/makefile.in

@@ -55,7 +55,7 @@ define cd_test
 
 endef
 
-test:
+test: fitpot
 	@$(foreach x,$(wildcard $(fitpot_examples)),$(call cd_test, ${x}))
 
 

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
 

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.

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

nappy/nn/make_params_NN.py

@@ -18,15 +18,16 @@
 import json
 
 _descfname = 'in.params.desc'
-_paramfname = 'in.params.NN2'
+_paramfname = 'in.params.'
 _logfname = 'log.make_params_NN'
 
 #...initial range of parameters
 _pmin = -20.0
 _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
 

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,...