Merge pull request #457 from abacusmodeling/develop

merge scalapack_gvx bugfix and related DFTU-bugfix

Author: dyzheng

source/input.cpp

@@ -1997,9 +1997,9 @@ bool Input::Read(const std::string &fn)
 			exit(0);
 		}
 
-		if(strcmp("genelpa", ks_solver.c_str())!=0)
+		if(strcmp("genelpa", ks_solver.c_str())!=0 && strcmp(ks_solver.c_str(),"scalapack_gvx")!=0 )
 		{
-			std::cout << " WRONG ARGUMENTS OF ks_solver in DFT+U routine, only genelpa is support " << std::endl;
+			std::cout << " WRONG ARGUMENTS OF ks_solver in DFT+U routine, only genelpa and scalapack_gvx are supportted " << std::endl;
 			exit(0);
 		}
 

source/src_lcao/DM_gamma.cpp

@@ -267,10 +267,7 @@ void Local_Orbital_Charge::allocate_gamma(const Grid_Technique &gt)
     setAlltoallvParameter(GlobalC::ParaO.comm_2D, GlobalC::ParaO.blacs_ctxt, GlobalC::ParaO.nb);
 
 	// Peize Lin test 2019-01-16
-    if (GlobalV::KS_SOLVER=="genelpa") //LiuXh add 2021-09-06, clear memory, _2d only used in genelpa solver
-    {
-        wfc_dm_2d.init();
-    }
+    wfc_dm_2d.init();
 
 	if(GlobalC::wf.start_wfc=="file")
 	{

tools/opt_orb_pytorch_dpsi/IO/cal_weight.py

@@ -0,0 +1,70 @@
+import IO.read_istate
+import torch
+import re
+import functools
+import operator
+
+def cal_weight(info_weight, flag_same_band, stru_file_list=None):
+	""" weight[ist][ib] """
+
+	if "bands_file" in info_weight.keys():
+		if "bands_range" in info_weight.keys():
+			raise IOError('"bands_file" and "bands_range" only once')
+
+		weight = []														# weight[ist][ib]
+		for weight_stru, file_name in zip(info_weight["stru"], info_weight["bands_file"]):
+			occ = IO.read_istate.read_istate(file_name)
+			weight += [occ_k * weight_stru for occ_k in occ]
+
+	elif "bands_range" in info_weight.keys():
+		k_weight = read_k_weight(stru_file_list)						# k_weight[ist][ik]
+		nbands = read_nbands(stru_file_list)							# nbands[ist]
+
+		st_weight = []													# st_weight[ist][ib]
+		for weight_stru, bands_range, nbands_ist in zip(info_weight["stru"], info_weight["bands_range"], nbands):
+			st_weight_tmp = torch.zeros((nbands_ist,))
+			st_weight_tmp[:bands_range] = weight_stru
+			st_weight.append( st_weight_tmp )
+
+		weight = []														# weight[ist][ib]
+		for ist,_ in enumerate(k_weight):
+			for ik,_ in enumerate(k_weight[ist]):
+				weight.append(st_weight[ist] * k_weight[ist][ik])
+
+	else:
+		raise IOError('"bands_file" and "bands_range" must once')
+
+
+	if not flag_same_band:
+		for ist,_ in enumerate(weight):
+			weight[ist] = torch.tensordot(weight[ist], weight[ist], dims=0)
+
+
+	normalization = functools.reduce(operator.add, map(torch.sum, weight), 0)
+	weight = list(map(lambda x:x/normalization, weight))
+
+	return weight
+
+
+def read_k_weight(stru_file_list):
+	""" weight[ist][ik] """
+	weight = []												# weight[ist][ik]
+	for file_name in stru_file_list:
+		weight_k = []										# weight_k[ik]
+		with open(file_name,"r") as file:
+			data = re.compile(r"<WEIGHT_OF_KPOINTS>(.+)</WEIGHT_OF_KPOINTS>", re.S).search(file.read()).group(1).split("\n")
+			for line in data:
+				line = line.strip()
+				if line:
+					weight_k.append(float(line.split()[-1]))
+		weight.append(weight_k)
+	return weight
+
+
+def read_nbands(stru_file_list):
+	""" nbands[ib] """
+	nbands = []
+	for file_name in stru_file_list:
+		with open(file_name,"r") as file:
+			nbands.append(int(re.compile(r"(\d+)\s+nbands").search(file.read()).group(1)))
+	return nbands

tools/opt_orb_pytorch_dpsi/IO/func_C.py

@@ -0,0 +1,87 @@
+from util import *
+import torch
+import numpy as np
+
+def random_C_init(info):
+	""" C[it][il][ie,iu]	<jY|\phi> """
+	C = dict()
+	for it in info.Nt_all:
+		C[it] = ND_list(info.Nl[it])
+		for il in range(info.Nl[it]):
+			C[it][il] = torch.tensor(np.random.uniform(-1,1, (info.Ne[it], info.Nu[it][il])), dtype=torch.float64, requires_grad=True)
+	return C
+	
+	
+	
+def read_C_init(file_name,info):
+	""" C[it][il][ie,iu]	<jY|\phi> """
+	C = random_C_init(info)
+
+	with open(file_name,"r") as file:
+	
+		for line in file:
+			if line.strip() == "<Coefficient>":	
+				line=None
+				break
+		ignore_line(file,1)
+	
+		C_read_index = set()
+		while True:
+			line = file.readline().strip()
+			if line.startswith("Type"):
+				it,il,iu = list(map(int,file.readline().split()));	
+				it=info.Nt_all[it-1];	iu-=1
+				C_read_index.add((it,il,iu))
+				line = file.readline().split()
+				for ie in range(info.Ne[it]):
+					if not line:	line = file.readline().split()
+					C[it][il].data[ie,iu] = float(line.pop(0))
+			elif line.startswith("</Coefficient>"):
+				break;
+			else:
+				raise IOError("unknown line in read_C_init "+file_name+"\n"+line)
+	return C, C_read_index
+
+	
+	
+def copy_C(C,info):
+	C_copy = dict()
+	for it in info.Nt_all:
+		C_copy[it] = ND_list(info.Nl[it])
+		for il in range(info.Nl[it]):
+			C_copy[it][il] = C[it][il].clone()
+	return C_copy
+	
+	
+	
+def write_C(file_name,info,C):
+	with open(file_name,"w") as file:
+		print("<Coefficient>", file=file)
+		print("\tTotal number of radial orbitals.", file=file)
+		for it,C_t in C.items():
+			for il,C_tl in enumerate(C_t):
+				for iu in range(C_tl.size()[1]):
+					print("\tType\tL\tZeta-Orbital", file=file)
+					print(f"\t  {info.Nt_all.index(it)+1} \t{il}\t    {iu+1}", file=file)
+					for ie in range(C_tl.size()[0]):
+						print("\t", C_tl[ie,iu].item(), file=file)
+		print("</Coefficient>", file=file)
+	
+	
+#def init_C(info):
+#	""" C[it][il][ie,iu] """
+#	C = ND_list(max(info.Nt))
+#	for it in range(len(C)):
+#		C[it] = ND_list(info.Nl[it])
+#		for il in range(info.Nl[it]):
+#			C[it][il] = torch.autograd.Variable( torch.Tensor( info.Ne, info.Nu[it][il] ), requires_grad = True )
+#			
+#	with open("C_init.dat","r") as file:
+#		line = []
+#		for it in range(len(C)):
+#			for il in range(info.Nl[it]):
+#				for i_n in range(info.Nu[it][il]):
+#					for ie in range(info.Ne[it]):
+#						if not line:	line=file.readline().split()
+#						C[it][il].data[ie,i_n] = float(line.pop(0))
+#	return C

tools/opt_orb_pytorch_dpsi/IO/print_QSV.py

@@ -0,0 +1,29 @@
+def print_V(V,file_name):
+	""" V[ist][ib] """
+	with open(file_name,"w") as file:
+		for V_s in V:
+			for V_sb in V_s:
+				print(1-V_sb.item(),end="\t",file=file)
+			print(file=file)
+			
+def print_S(S,file_name):
+	""" S[ist][it1,it2][il1][il2][ia1*im1*in1,ia2*im2*in2] """
+	with open(file_name,"w") as file:
+		for ist,S_s in enumerate(S):
+			for (it1,it2),S_tt in S_s.items():
+				for il1,S_ttl in enumerate(S_tt):
+					for il2,S_ttll in enumerate(S_ttl):
+						print(ist,it1,it2,il1,il2,file=file)
+						print(S_ttll.real.numpy(),file=file)
+						print(S_ttll.imag.numpy(),"\n",file=file)
+						
+def print_Q(Q,file_name):
+	""" Q[ist][it][il][ib,ia*im*iu] """
+	with open(file_name,"w") as file:
+		for ist,Q_s in enumerate(Q):
+			for it,Q_st in Q_s.items():
+				for il,Q_stl in enumerate(Q_st):
+					print(ist,it,il,file=file)
+					print(Q_stl.real.numpy(),file=file)
+					print(Q_stl.imag.numpy(),"\n",file=file)
+					

tools/opt_orb_pytorch_dpsi/IO/print_orbital.py

@@ -0,0 +1,43 @@
+def print_orbital(orb,info):
+	""" orb[it][il][iu][r] """
+	for it,orb_t in orb.items():
+		with open("orb_{0}.dat".format(it),"w") as file:
+			print_orbital_head(file,info,it)
+			for il,orb_tl in enumerate(orb_t):
+				for iu,orb_tlu in enumerate(orb_tl):
+					print("""                Type                   L                   N""",file=file)
+					print("""                   0                   {0}                   {1}""".format(il,iu),file=file)
+					for ir,orb_tlur in enumerate(orb_tlu):
+						print(orb_tlur,end="\t",file=file)
+						if ir%4==3:	print(file=file)
+					print(file=file)
+					
+					
+def plot_orbital(orb,Rcut,dr):
+	for it,orb_t in orb.items():
+		with open("orb_{0}_plot.dat".format(it),"w") as file:
+			Nr = int(Rcut[it]/dr[it])+1
+			for ir in range(Nr):
+				print(ir*dr[it],end="\t",file=file)
+				for il,orb_tl in enumerate(orb_t):
+					for orb_tlu in orb_tl:
+						print(orb_tlu[ir],end="\t",file=file)
+				print(file=file)
+				
+				
+def print_orbital_head(file,info,it):
+	print( "---------------------------------------------------------------------------", file=file )
+	print( "Element                     {0}".format(it), file=file )
+	print( "Energy Cutoff(Ry)           {0}".format(info.Ecut[it]), file=file )
+	print( "Radius Cutoff(a.u.)         {0}".format(info.Rcut[it]), file=file )
+	print( "Lmax                        {0}".format(info.Nl[it]-1), file=file )
+	l_name = ["S","P","D"]+list(map(chr,range(ord('F'),ord('Z')+1)))
+	for il,iu in enumerate(info.Nu[it]):
+		print( "Number of {0}orbital-->       {1}".format(l_name[il],iu), file=file )
+	print( "---------------------------------------------------------------------------", file=file )
+	print( "SUMMARY  END", file=file )
+	print( file=file )
+	print( "Mesh                        {0}".format(int(info.Rcut[it]/info.dr[it])+1), file=file )
+	print( "dr                          {0}".format(info.dr[it]), file=file )
+
+