关于 energy_grand 的计算问题

[IMPLICIT-NONE]

您好，我目前在学习 GC-DFT，还是新手，看了您的代码 vaspFCP.py，有些不明白的地方。在 vaspFCP.py 中，energy_grand 是这么计算的：energy_grand1=float(line.split()[-2])+(self.U+self.workSHE+self.fermishift)*(self.Nelect-self.Nelect0)。我的疑问是：

1. self.U 是设定的目标电势，这里是不是应该用当前计算下的电势值（Ucal）？
2. 如果是用 Ucal 的话，Ucal=-(self.fermi + self.fermishift)-self.workSHE，带入 energy_grand1 的计算式子中，self.U+self.workSHE+self.fermishift = self.fermi，也就是说 energy_grand1 = =float(line.split()[-2])+ self.fermi * (self.Nelect-self.Nelect0)。然而，这里用的 self.fermi 是 vasp 给出的费米能级，是还没有调整为 SHE 标准下的数值，为此我感到疑惑。这里是不是应该使用 SHE 标准下的费米能级？即 energy_grand1=float(line.split()[-2])+(self.fermi + self.fermishift)*(self.Nelect-self.Nelect0) 。

P.S.：代码的输出有点乱，建议考虑格式化输出。这是我自己改的部分，供参考：

1. 230 行：
   print("%10s %12s %12s %12s %12s %12s %12s %12s" % ('loop', 'Nelect', 'feimi', 'fermishift', 'energyclose', 'energygrand', 'Ucal', 'C'), file = f)
   #f.write('loop'+' '+'Nelect' +' '+'feimi'+' '+'fermishift'+ ' ' +'energyclose' +' '+'energygrand'+' '+'Ucal'+' '+'C'+'\n')
2. 453 行：
   print("%10s %12s %12s %12s %12s %12s %12s %12s %12s" % ('loop', 'NELECT', 'fermi', 'fermishift', 'U', 'conv', 'G', 'G_grand', 'time'), file = f)
   #f.write('loop'+' '+'NELECT'+' ' +'fermi'+' '+'fermishift'+' ' + 'U' + ' ' + 'conv' +' '+'G'+' '+'G_grand'+' '+ 'time'+'\n')
3. 484 行：
   print("%10.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f" % (self.FCPloop, self.Nelect, self.fermi, self.fermishift, Ucal, conv, self.energy_close, self.energy_free, endcal-startcal), file = f)
   #f.write(str(self.FCPloop)+' '+str(self.Nelect) +' '+str(self.fermi)+' '+str(self.fermishift)+' '+ str(Ucal)+' ' + str(conv)+ ' ' +str(self.energy_close) +' '+str(self.energy_free)+' '+str(endcal-startcal)+'\n')
4. 572 行：
   print("%10.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f" % (self.FCPloop, self.Nelect, self.fermi, self.fermishift, Ucal, conv, self.energy_close, self.energy_free, endcal-startcal), file = f)
   #f.write(str(self.FCPloop)+' '+str(self.Nelect) +' '+str(self.fermi)+' '+str(self.fermishift)+' '+ str(Ucal)+' ' + str(conv)+ ' ' +str(self.energy_close) +' '+str(self.energy_free)+' '+str(endcal-startcal)+'\n')
5. 604 行：
   print("%10.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f" % (self.FCPloop, Nelectold, self.fermi, self.fermishift, self.energy_close, self.energy_free, Ucal, self.C), file = f)
   #f.write(str(self.FCPloop)+' '+str(Nelectold) +' '+str(self.fermi)+' '+str(self.fermishift)+ ' ' +str(self.energy_close) +' '+str(self.energy_free)+' '+str(Ucal)+' '+str(self.C)+'\n')

[hellozhaoming]

你好，由于我在做其他项目，暂时没有太多时间去维护。感谢你提出的格式化输出的建议。

1.巨正则系综中，粒子可以在系统和环境之间发生交换。巨正则能量的定义中化学势是环境的化学势，这个化学势相当于能量参考点，所以应该用self.U，而不是Ucal。如果用Ucal，那么能量参考点就在变化，能量的对比就没有意义。

2. self.U+self.workSHE+self.fermishift 这个计算式中理论上不应该包含self.fermishift这一项。但是vasp计算能量时电子的能量参考点是整个晶胞的平均静电势而不是真空静电势。当体系为电中性时没有影响，但是当体系带有净电荷时，真空静电势会带来误差，即float(line.split()[-2])这一项需要对真空能级进行矫正，矫正系数为self.fermishift*(self.Nelect-self.Nelect0)。矫正过后能量参考点均为真空。self.U+self.workSHE的能量参考点也是真空。
   当然，也可以用其他能量参考点定义巨正则能量。比如用SHE作为参考点。float(line.split()[-2])这一项的电子能也需要矫正到SHE参考点下的能量，所以还是会包含self.workSHE这一项。

如果你不确定哪个公式是正确的，你可以改变体系的电荷数，分别用每种能量公式计算电荷与能量的关系，观察能量极小值点对应的电荷是否是目标电势对应的电荷。

[IMPLICIT-NONE]

非常感谢您的解答。

我之前应该是把 energy_grand 和恒电势模型中的吉布斯自由能搞混了。 energy_grand 是设定 self.U 下的极小值，用来求解 self.U 的，吉布斯自由能是用来描绘催化反应的反应路径的。

我还有一些问题想请教一下：

1. 对于 FCP cycle 中不容易收敛的体系，使用 adaptive_lr = True 似乎能有效改善收敛性。为了改善 FCP 收敛性，对于大多数体系，在计算的时候是否可以直接开启 adaptive_lr = True？
2. 我在尝试除了 Newton-fitting 以外的 FCPmethod 的时候遇到了报错，看了 vaspFCP.py，似乎其他的方法还没有补充进来？

另外我还有一些建议，供您参考：

1. 当使用 PREC = N 的时候，在个别体系的自洽计算中，可能会遇到 FERMI_SHIFT 全部等于 0 的情况，导致 FCP 计算无法获取 E_fs 而报错退出，计算中断。FERMI_SHIFT = 0 可能是 vaspsol 的代码的一个 bug。使用 PREC = A （其实就是增加 NGX，NGY，NGZ）可以避免 FERMI_SHIFT = 0 的情况出现。可以考虑在 readme 文件里面建议使用 PREC = A。
2. 磁性体系的计算经常会遇到计算不收敛或者是收敛到错误的磁态上面，这可能会影响 FCP 计算的结果。因此，对于磁性体系，检查计算结果（能量以及磁矩）是很有必要的。目前的 vaspFCP.py 会覆盖 vasp.out 和 tmp-log-FCP.txt 文件，会给检查带来不便。可以考虑把所有的 vasp 输出信息以及 tmp-log-FCP.txt 信息都保留下来。这是我自己用的代码（daoxiong add start 和 daoxiong add end标记），供参考：

               while self.FCPloop < loopmax:
                startcal=time.time()
                self.FCPloop += 1
                self.write_input(self.atoms, properties, system_changes)
                if self.Nelect != None:
                    os.system('echo NELECT='+ str(self.Nelect) + ' >> ' + self.directory +'/INCAR')
                with self._txt_outstream() as out:
                    errorcode = self._run(command=command,
                                          out=out,
                                          directory=self.directory)

                if errorcode:
                    raise calculator.CalculationFailed(
                        '{} in {} returned an error: {:d}'.format(
                            self.name, self.directory, errorcode))

                #### daoxiong add start ####
                # Post-processing and read results from calculation
                os.system('echo -------- FCPloop = ' + str(self.FCPloop) + ' -------- >> ' + self.directory + '/log.vasp-fcp')
                os.system('cat  ' + self.directory + '/' + self.txt + ' >> ' + self.directory + '/log.vasp-fcp')
                #### daoxiong add end ####

                self.update_atoms(atoms)
                self.read_results()
                Ucal=-(self.fermi + self.fermishift)-self.workSHE
                conv=self.U + self.workSHE + self.fermi + self.fermishift
                endcal=time.time()
                with open(LogPath+ '/tmp-log-FCP.txt', mode='a',encoding='utf-8') as f:
                    print("%4.0f  %11.6f  %11.6f  %11.6f  %7.3f  %11.6f  %11.6f  %11.6f  %11.6f  %11.6f  %7.0f" % (self.FCPloop, self.Nelect, self.fermi, self.fermishift, self.workSHE, Ucal, conv, self.energy_close, self.energy_free, self.C, endcal-startcal), file = f)
                    #f.write(str(self.FCPloop)+' '+str(self.Nelect) +' '+str(self.fermi)+' '+str(self.fermishift)+' '+ str(Ucal)+' ' + str(conv)+ ' ' +str(self.energy_close) +' '+str(self.energy_free)+' '+str(endcal-startcal)+'\n')

                Ctemp= Cevalue(self.Nelect)
                if Ctemp > 0:
                    self.C=Ctemp
                if convold != None:
                    if convold*conv < 0:
                        lrcount+=1
                convold = conv
                Nelectold = self.Nelect

                if self.FCPloop== loopmax:
                    break

                if self.FCPmethod == 'Newton-fitting':
                    if self.adaptive_lr==True:
                        if lrcount >1:
                            lr=lr-0.5**lrcount
                    #print(str(Nelectold),str(lr),str(convold),str(self.C))
                    self.Nelect=Nelectold-lr*convold*self.C  #lr is learning rate

                if abs(self.Nelect-Nelectold)/self.Nelect > 0.05:
                    os.system('rm '+self.directory +'/WAVECAR')

                if abs(conv)<FCPconv:
                    break

        with open('./log-' + str(1) + '.txt', mode='a',encoding='utf-8') as f:
            print("%7.6f  %11.6f  %11.6f  %11.6f  %11.6f  %11.6f  %11.6f  %11.6f  %11.6f" % (self.FCPloop, Nelectold, self.fermi, self.fermishift, self.workSHE, self.energy_close, self.energy_free, Ucal, self.C), file = f)
            #f.write(str(self.FCPloop)+' '+str(Nelectold) +' '+str(self.fermi)+' '+str(self.fermishift)+ ' ' +str(self.energy_close) +' '+str(self.energy_free)+' '+str(Ucal)+' '+str(self.C)+'\n')

        #### daoxiong add start ####
        os.system('cat  ' + self.directory + '/tmp-log-FCP.txt  >> ' + self.directory + '/log.tmp-log-FCP')
        #### daoxiong add end ####

[hellozhaoming]

adaptive_lr = True对于一些比较好收敛的体系会严重拖慢速度，所以默认是关闭的。建议遇到不收敛时去看一下log文件，作Nelect和E_grand的图。对于磁性体系，电子态的能量极小值点不唯一，可能在不同电子态之间来回跳跃导致不收敛，即使adaptive_lr = True使得体系收敛了，也很难确保这个收敛的态就是电子态的基态，建议仔细检查。

Newton-fitting 以外的方法收敛速度很慢，我在新版本中去掉了

FERMI_SHIFT 全部等于 0 有很多种可能的原因，比如真空层厚度不够，电子收敛精度不够都有可能发生。

非常感谢你的建议，我迭代了一个新版本。这个版本应该不止可以应用于vasp，还可以自定义其他first princple的calculator，只要自己写get_fermi()和get_fermi_shift()函数。当然我在忙其他项目，暂时没时间维护。如果你有空的话可以新建一个github分支加进去