Python numpy.savetxt() Examples

def test_bse_gto_vs_nao_inter_0082(self):
    """ Interacting case """
    #dm1 = gto_mf.make_rdm1()
    #o1 = gto_mf.get_ovlp()
    #print(__name__, 'dm1*o1', (dm1*o1).sum())
    nao_td = bse_iter(mf=gto_mf, gto=mol, verbosity=0, xc_code='GW', perform_gw=True)
    
    #dm2 = nao_td.make_rdm1()
    #o2 = nao_td.get_ovlp()
    #n = nao_td.norbs
    #print(__name__, 'dm2*o2', (dm2.reshape((n,n))*o2).sum())
    
    omegas = np.linspace(0.0,2.0,450)+1j*0.04
    p_iter = -nao_td.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*27.2114, p_iter])
    np.savetxt('be.bse_iter.omega.inter.ave.txt', data.T, fmt=['%f','%f']) 

def test_0150_bse_h2o_uhf_rpa(self):
    """ Interacting case """
    mol=gto.M(verbose=0,atom='O 0 0 0;H 0 0.489 1.074;H 0 0.489 -1.074',basis='cc-pvdz',)
    gto_mf = scf.UKS(mol)
    gto_mf.xc = 'hf'
    gto_mf.kernel()
    gto_td = tddft.dRPA(gto_mf)
    gto_td.nstates = 95
    gto_td.kernel()

    omegas = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    p_ave = -polariz_freq_osc_strength(gto_td.e, gto_td.oscillator_strength(), omegas).imag
    data = np.array([omegas.real*HARTREE2EV, p_ave])
    np.savetxt('test_0150_bse_h2o_uhf_rpa_pyscf.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0150_bse_h2o_uhf_rpa_pyscf.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))
    
    nao_td  = bse_iter(mf=gto_mf, gto=mol, verbosity=0, xc_code='RPA')

    polariz = -nao_td.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*HARTREE2EV, polariz])
    np.savetxt('test_0150_bse_h2o_uhf_rpa_nao.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0150_bse_h2o_uhf_rpa_nao.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3)) 

def test_151_bse_h2b_uhf_rpa(self):
    """ This  """
    mol = gto.M(verbose=1,atom='B 0 0 0; H 0 0.489 1.074; H 0 0.489 -1.074',basis='cc-pvdz',spin=3)

    gto_mf = scf.UKS(mol)
    gto_mf.xc = 'hf'
    gto_mf.kernel()
    gto_td = tddft.dRPA(gto_mf)
    gto_td.nstates = 190
    gto_td.kernel()

    omegas = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    p_ave = -polariz_freq_osc_strength(gto_td.e, gto_td.oscillator_strength(), omegas).imag
    data = np.array([omegas.real*HARTREE2EV, p_ave])
    np.savetxt('test_0151_bse_h2b_uhf_rpa_pyscf.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0151_bse_h2b_uhf_rpa_pyscf.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))
    
    nao_td  = bse_iter(mf=gto_mf, gto=mol, verbosity=0, xc_code='RPA')

    polariz = -nao_td.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*HARTREE2EV, polariz])
    np.savetxt('test_0151_bse_h2b_uhf_rpa_nao.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0151_bse_h2b_uhf_rpa_nao.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3)) 

def test_qchem_irf(self):
    """ Test """
    gto_hf = dft.RKS(mol)
    gto_hf.kernel()
    nocc = mol.nelectron//2
    nmo = gto_hf.mo_energy.size
    nvir = nmo-nocc    
    gto_td = tddft.dRPA(gto_hf)
    gto_td.nstates = min(100, nocc*nvir)
    gto_td.kernel()
    
    gto_gw = GW(gto_hf, gto_td)
    gto_gw.kernel()
    
    nao_td  = tddft_iter(mf=gto_hf, gto=mol, xc_code='RPA')
    eps = 0.02
    omegas = np.arange(0.0,2.0,eps/2.0)+1j*eps
    p_iter = -nao_td.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*27.2114, p_iter])
    np.savetxt('NH.tddft_iter_rpa.omega.inter.pav.txt', data.T, fmt=['%f','%f'])
    
    np.set_printoptions(linewidth=180)
    #qrf = qchem_inter_rf(mf=gto_hf, gto=mol, pb_algorithm='fp', verbosity=1) 

def pred_test(testing_data, exe, param_list=None, save_path=""):
    ret = numpy.zeros((testing_data.shape[0], 2))
    if param_list is None:
        for i in range(testing_data.shape[0]):
            exe.arg_dict['data'][:] = testing_data[i, 0]
            exe.forward(is_train=False)
            ret[i, 0] = exe.outputs[0].asnumpy()
            ret[i, 1] = numpy.exp(exe.outputs[1].asnumpy())
        numpy.savetxt(save_path, ret)
    else:
        for i in range(testing_data.shape[0]):
            pred = numpy.zeros((len(param_list),))
            for j in range(len(param_list)):
                exe.copy_params_from(param_list[j])
                exe.arg_dict['data'][:] = testing_data[i, 0]
                exe.forward(is_train=False)
                pred[j] = exe.outputs[0].asnumpy()
            ret[i, 0] = pred.mean()
            ret[i, 1] = pred.std()**2
        numpy.savetxt(save_path, ret)
    mse = numpy.square(ret[:, 0] - testing_data[:, 0] **3).mean()
    return mse, ret 

def test_sv_after_gpaw(self):
    """ init ao_log with radial orbitals from GPAW """
    if calc is None: return

    self.assertTrue(hasattr(calc, 'setups'))
    sv = mf(gpaw=calc, gen_pb=False)
    self.assertEqual(sv.ao_log.nr, 1024)
    over = sv.overlap_coo().toarray()
    error = sv.hsx.check_overlaps(over)
    self.assertLess(error, 1e-4)
    #print("overlap error: ", error/over.size)

    #print("Pyscf    gpaw")
    #for py, gp in zip(over[:, 0], sv.wfsx.overlaps[:, 0]):
    #    print("{0:.5f}    {1:.5f}".format(py, gp))
    #np.savetxt("Pyscf.overlaps", over)
    #np.savetxt("gpaw.overlaps", sv.wfsx.overlaps,) 

def write_param(self, xml_filename='gap.xml'):
        """
        Write xml file to perform lammps calculation.

        Args:
            xml_filename (str): Filename to store xml formatted parameters.
        """
        if not self.param:
            raise RuntimeError("The xml and parameters should be provided.")
        tree = self.param.get('xml')
        root = tree.getroot()
        gpcoordinates = list(root.iter('gpCoordinates'))[0]
        param_filename = "{}.soapparam".format(self.name)
        gpcoordinates.set('sparseX_filename', param_filename)
        np.savetxt(param_filename, self.param.get('param'), fmt='%.20e')
        tree.write(xml_filename)
        pair_coeff = self.pair_coeff.format(xml_filename,
                                            '\"Potential xml_label={}\"'.
                                            format(self.param.get('potential_label')),
                                            self.specie.Z)
        ff_settings = [self.pair_style, pair_coeff]
        return ff_settings 

def print_mutation(hyp, results, bucket=''):
    # Print mutation results to evolve.txt (for use with train.py --evolve)
    a = '%10s' * len(hyp) % tuple(hyp.keys())  # hyperparam keys
    b = '%10.3g' * len(hyp) % tuple(hyp.values())  # hyperparam values
    c = '%10.3g' * len(results) % results  # results (P, R, mAP, F1, test_loss)
    print('\n%s\n%s\nEvolved fitness: %s\n' % (a, b, c))

    if bucket:
        os.system('gsutil cp gs://%s/evolve.txt .' % bucket)  # download evolve.txt

    with open('evolve.txt', 'a') as f:  # append result
        f.write(c + b + '\n')
    x = np.unique(np.loadtxt('evolve.txt', ndmin=2), axis=0)  # load unique rows
    np.savetxt('evolve.txt', x[np.argsort(-fitness(x))], '%10.3g')  # save sort by fitness

    if bucket:
        os.system('gsutil cp evolve.txt gs://%s' % bucket)  # upload evolve.txt 

def run(self):
        """
        resize images then write manifest files to disk.
        """
        self.write_label()
        self.collectdata()

        records = [(fname, tgt)
                   for fname, tgt in self.trainpairlist.items()]
        np.savetxt(self.manifests['train'], records, fmt='%s,%s')

        records = [(fname, tgt)
                   for fname, tgt in self.valpairlist.items()]
        np.savetxt(self.manifests['val'], records, fmt='%s,%s')

        records = [(fname, tgt)
                   for fname, tgt in self.testpairlist.items()]
        np.savetxt(self.manifests['test'], records, fmt='%s,%s') 

def extract_feature(model, dataloader, save_path, load_from_disk=True, model_path=''):
    if load_from_disk:
        model = models.Network(base_net=args.model_name,
                               n_class=args.num_class)
        model.load_state_dict(torch.load(model_path))
        model = model.to(DEVICE)
    model.eval()
    correct = 0
    fea_all = torch.zeros(1,1+model.base_network.output_num()).to(DEVICE)
    with torch.no_grad():
        for inputs, labels in dataloader:
            inputs, labels = inputs.to(DEVICE), labels.to(DEVICE)
            feas = model.get_features(inputs)
            labels = labels.view(labels.size(0), 1).float()
            x = torch.cat((feas, labels), dim=1)
            fea_all = torch.cat((fea_all, x), dim=0)
            outputs = model(inputs)
            preds = torch.max(outputs, 1)[1]
            correct += torch.sum(preds == labels.data.long())
        test_acc = correct.double() / len(dataloader.dataset)
    fea_numpy = fea_all.cpu().numpy()
    np.savetxt(save_path, fea_numpy[1:], fmt='%.6f', delimiter=',')
    print('Test acc: %f' % test_acc)

# You may want to classify with 1nn after getting features 

def test_139_h2o_uks_nonin_pb(self):
    """ This  """
    mol = gto.M(verbose=1,atom='O 0 0 0; H 0 0.489 1.074; H 0 0.489 -1.074',basis='cc-pvdz')
    gto_mf = dft.UKS(mol)
    gto_mf.kernel()
    nao_mf = bse_iter(gto=mol, mf=gto_mf)
    comega = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    pnonin = -nao_mf.polariz_nonin_ave_matelem(comega).imag
    data = np.array([comega.real*HARTREE2EV, pnonin])
    np.savetxt('test_139_h2o_uks_nonin_matelem.txt', data.T, fmt=['%f','%f'])
    #data_ref = np.loadtxt('test_139_h2o_uks_nonin_matelem.txt-ref').T
    #self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))

    pnonin = -nao_mf.comp_polariz_nonin_ave(comega).imag
    data = np.array([comega.real*HARTREE2EV, pnonin])
    np.savetxt('test_139_h2o_uks_nonin_nao.txt', data.T, fmt=['%f','%f'])
    #data_ref = np.loadtxt('test_139_h2o_uks_nonin_nao.txt-ref').T
    #self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3)) 

def test_152_bse_h2b_uks_pz(self):
    """ This  """
    mol = gto.M(verbose=1,atom='B 0 0 0; H 0 0.489 1.074; H 0 0.489 -1.074',basis='cc-pvdz',spin=3)

    gto_mf = scf.UKS(mol)
    gto_mf.kernel()
    gto_td = tddft.TDDFT(gto_mf)
    gto_td.nstates = 190
    gto_td.kernel()

    omegas = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    p_ave = -polariz_freq_osc_strength(gto_td.e, gto_td.oscillator_strength(), omegas).imag
    data = np.array([omegas.real*HARTREE2EV, p_ave])
    np.savetxt('test_0152_bse_h2b_uks_pz_pyscf.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0152_bse_h2b_uks_pz_pyscf.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))
    
    nao_td  = bse_iter(mf=gto_mf, gto=mol, verbosity=0)

    polariz = -nao_td.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*HARTREE2EV, polariz])
    np.savetxt('test_0152_bse_h2b_uks_pz_nao.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0152_bse_h2b_uks_pz_nao.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3)) 

def test_bse_iter_nonin(self):
    """ Compute polarization with LDA TDDFT  """
    from timeit import default_timer as timer
    
    dname = os.path.dirname(os.path.abspath(__file__))
    bse = bse_iter(label='water', cd=dname, iter_broadening=1e-2, xc_code='RPA', verbosity=0)
    omegas = np.linspace(0.0,2.0,500)+1j*bse.eps
    
    pxx = np.zeros(len(omegas))
    for iw,omega in enumerate(omegas):
      for ixyz in range(1):
        vab = bse.apply_l0(bse.dip_ab[ixyz], omega)
        pxx[iw] = pxx[iw] - (vab.imag*bse.dip_ab[ixyz].reshape(-1)).sum()
        
    data = np.array([omegas.real*27.2114, pxx])
    np.savetxt('water.bse_iter_rpa.omega.nonin.pxx.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt(dname+'/water.bse_iter_rpa.omega.nonin.pxx.txt-ref')
    #print('    bse.l0_ncalls ', bse.l0_ncalls)
    self.assertTrue(np.allclose(data_ref,data.T, rtol=1.0, atol=1e-05)) 

def test_0146_bse_h2o_rhf_rpa(self):
    """ Interacting case """
    mol=gto.M(verbose=0,atom='O 0 0 0;H 0 0.489 1.074;H 0 0.489 -1.074',basis='cc-pvdz',)
    gto_hf = scf.RKS(mol)
    gto_hf.xc = 'hf'
    gto_hf.kernel()
    gto_td = tddft.dRPA(gto_hf)
    gto_td.nstates = 95
    gto_td.kernel()

    omegas = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    #p_ave = -polariz_inter_ave(gto_hf, mol, gto_td, omegas).imag
    p_ave = -polariz_freq_osc_strength(gto_td.e, gto_td.oscillator_strength(), omegas).imag
    data = np.array([omegas.real*HARTREE2EV, p_ave])
    np.savetxt('test_0146_bse_h2o_rhf_rpa_pyscf.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0146_bse_h2o_rhf_rpa_pyscf.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))
    
    nao_td  = bse_iter(mf=gto_hf, gto=mol, verbosity=0, xc_code='RPA')

    p_iter = -nao_td.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*HARTREE2EV, p_iter])
    np.savetxt('test_0146_bse_h2o_rhf_rpa_nao.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0146_bse_h2o_rhf_rpa_nao.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3)) 

def test_gw(self):
    """ This is GW """
    mol = gto.M( verbose = 1, atom = '''H 0 0 0;  H 0.17 0.7 0.587''', basis = 'cc-pvdz',)
    #mol = gto.M( verbose = 0, atom = '''H 0.0 0.0 -0.3707;  H 0.0 0.0 0.3707''', basis = 'cc-pvdz',)
    gto_mf = scf.RHF(mol)
    gto_mf.kernel()
    #print('gto_mf.mo_energy:', gto_mf.mo_energy)
    b = bse_iter(mf=gto_mf, gto=mol, perform_gw=True, xc_code='GW', verbosity=0, nvrt=4)
    #self.assertAlmostEqual(b.mo_energy[0], -0.5967647)
    #self.assertAlmostEqual(b.mo_energy[1], 0.19072719)
    omegas = np.linspace(0.0,2.0,450)+1j*0.04
    p_iter = -b.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*27.2114, p_iter])
    np.savetxt('h2_gw_bse_iter.omega.inter.ave.txt', data.T)
    data_ref = np.loadtxt('h2_gw_bse_iter.omega.inter.ave.txt-ref').T
    #print(__name__, abs(data_ref-data).sum()/data.size)
    self.assertTrue(np.allclose(data_ref, data, 5))

    p_iter = -b.comp_polariz_nonin_ave(omegas).imag
    data = np.array([omegas.real*27.2114, p_iter])
    np.savetxt('h2_gw_bse_iter.omega.nonin.ave.txt', data.T) 

def test_0040_bse_rpa_nao(self):
    """ Compute polarization with RPA via 2-point non-local potentials (BSE solver)  """
    from timeit import default_timer as timer

    dname = os.path.dirname(os.path.abspath(__file__))
    bse = bse_iter(label='water', cd=dname, xc_code='RPA', verbosity=0)
    omegas = np.linspace(0.0,2.0,150)+1j*0.01
    #print(__name__, omegas.shape)
    
    pxx = np.zeros(len(omegas))
    for iw,omega in enumerate(omegas):
      for ixyz in range(1):
        dip = bse.dip_ab[ixyz]
        vab = bse.apply_l(dip, omega)
        pxx[iw] = pxx[iw] - (vab.imag*dip.reshape(-1)).sum()
        
    data = np.array([omegas.real*27.2114, pxx])
    np.savetxt('water.bse_iter_rpa.omega.inter.pxx.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt(dname+'/water.bse_iter_rpa.omega.inter.pxx.txt-ref')
    #print('    bse.l0_ncalls ', bse.l0_ncalls)
    self.assertTrue(np.allclose(data_ref,data.T, rtol=1.0, atol=1e-05)) 

def test_0162_bse_h2o_spin2_uhf_cis(self):
    """ This  """
    mol = gto.M(verbose=1,atom='O 0 0 0; H 0 0.489 1.074; H 0 0.489 -1.074',basis='cc-pvdz',spin=2)

    gto_mf = scf.UHF(mol)
    gto_mf.kernel()
    gto_td = tddft.TDDFT(gto_mf)
    gto_td.nstates = 190
    gto_td.kernel()

    omegas = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    p_ave = -polariz_freq_osc_strength(gto_td.e, gto_td.oscillator_strength(), omegas).imag
    data = np.array([omegas.real*HARTREE2EV, p_ave])
    np.savetxt('test_0162_bse_h2o_spin2_uhf_cis_pyscf.txt', data.T, fmt=['%f','%f'])
    #data_ref = np.loadtxt('test_0162_bse_h2o_spin2_uhf_cis_pyscf.txt-ref').T
    #self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))
    
    nao_td  = bse_iter(mf=gto_mf, gto=mol, verbosity=0, xc_code='CIS')

    polariz = -nao_td.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*HARTREE2EV, polariz])
    np.savetxt('test_0162_bse_h2o_spin2_uhf_cis_nao.txt', data.T, fmt=['%f','%f'])
    #data_ref = np.loadtxt('test_0162_bse_h2o_spin2_uhf_cis_nao.txt-ref').T
    #self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3), \
    #  msg="{}".format(abs(data_ref-data).sum()/data.size)) 

def test_x_zip_feature_na20_chain(self):
    """ This a test for compression of the eigenvectos at higher energies """
    dname = dirname(abspath(__file__))
    siesd = dname+'/sodium_20'
    x = td_c(label='siesta', cd=siesd,x_zip=True, x_zip_emax=0.25,x_zip_eps=0.05,jcutoff=7,xc_code='RPA',nr=128, fermi_energy=-0.0913346431431985)
    
    eps = 0.005
    ww = np.arange(0.0, 0.5, eps/2.0)+1j*eps
    data = np.array([ww.real*27.2114, -x.comp_polariz_inter_ave(ww).imag])
    fname = 'na20_chain.tddft_iter_rpa.omega.inter.ave.x_zip.txt'
    np.savetxt(fname, data.T, fmt=['%f','%f'])
    #print(__file__, fname)
    data_ref = np.loadtxt(dname+'/'+fname+'-ref')
    #print('    x.rf0_ncalls ', x.rf0_ncalls)
    #print(' x.matvec_ncalls ', x.matvec_ncalls)
    self.assertTrue(np.allclose(data_ref,data.T, rtol=1.0e-1, atol=1e-06)) 

def get_predict_labels():
    inputs = tf.placeholder("float", [None, 64, 64, 1])
    is_training = tf.placeholder("bool")
    prediction, _ = googlenet(inputs, is_training)
    predict_labels = tf.argmax(prediction, 1)
    sess = tf.Session()
    sess.run(tf.global_variables_initializer())
    saver = tf.train.Saver()
    data = sio.loadmat("../data/dataset.mat")
    testdata = data["test"] / 127.5 - 1.0
    testlabel = data["testlabels"]
    saver.restore(sess, "../save_para/.\\model.ckpt")
    nums_test = testlabel.shape[1]
    PREDICT_LABELS = np.zeros([nums_test])
    for i in range(nums_test // BATCH_SIZE):
        PREDICT_LABELS[i * BATCH_SIZE:i * BATCH_SIZE + BATCH_SIZE] = sess.run(predict_labels, feed_dict={inputs: testdata[i * BATCH_SIZE:i * BATCH_SIZE + BATCH_SIZE], is_training: False})
    PREDICT_LABELS[(nums_test // BATCH_SIZE - 1) * BATCH_SIZE + BATCH_SIZE:] = sess.run(predict_labels, feed_dict={inputs: testdata[(nums_test // BATCH_SIZE - 1) * BATCH_SIZE + BATCH_SIZE:], is_training: False})
    np.savetxt("../data/predict_labels.txt", PREDICT_LABELS) 

def test_161_bse_h2b_spin1_uhf_cis(self):
    """ This  """
    mol = gto.M(verbose=1,atom='B 0 0 0; H 0 0.489 1.074; H 0 0.489 -1.074',basis='cc-pvdz',spin=1)

    gto_mf = scf.UHF(mol)
    gto_mf.kernel()
    gto_td = tddft.TDHF(gto_mf)
    gto_td.nstates = 150
    gto_td.kernel()

    omegas = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    p_ave = -polariz_freq_osc_strength(gto_td.e, gto_td.oscillator_strength(), omegas).imag
    data = np.array([omegas.real*HARTREE2EV, p_ave])
    np.savetxt('test_0161_bse_h2b_spin1_uhf_cis_pyscf.txt', data.T, fmt=['%f','%f'])
    #data_ref = np.loadtxt('test_0159_bse_h2b_uhf_cis_pyscf.txt-ref').T
    #self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))
    
    nao_td  = bse_iter(mf=gto_mf, gto=mol, verbosity=0, xc_code='CIS')

    polariz = -nao_td.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*HARTREE2EV, polariz])
    np.savetxt('test_0161_bse_h2b_spin1_uhf_cis_nao.txt', data.T, fmt=['%f','%f'])
    #data_ref = np.loadtxt('test_0161_bse_h2b_spin1_uhf_cis_nao.txt-ref').T
    #self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3)) 

def generate_test_data():
    for str_problem in ["osy"]:
        problem = get_problem(str_problem)

        X = []

        # define a callback function that prints the X and F value of the best individual
        def my_callback(algorithm):
            pop = algorithm.pop
            _X = pop.get("X")[np.random.permutation(len(pop))[:10]]
            X.append(_X)

        minimize(problem,
                 method='nsga2',
                 method_args={'pop_size': 100},
                 termination=('n_gen', 100),
                 callback=my_callback,
                 pf=problem.pareto_front(),
                 disp=True,
                 seed=1)

        np.savetxt("%s.x" % str_problem, np.concatenate(X, axis=0), delimiter=",") 

def test_157_bse_h2o_rhf_cis(self):
    """ This  """
    mol = gto.M(verbose=1,atom='O 0 0 0; H 0 0.489 1.074; H 0 0.489 -1.074',basis='cc-pvdz')

    gto_mf = scf.RHF(mol)
    gto_mf.kernel()
    gto_td = tddft.TDDFT(gto_mf)
    gto_td.nstates = 190
    gto_td.kernel()

    omegas = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    p_ave = -polariz_freq_osc_strength(gto_td.e, gto_td.oscillator_strength(), omegas).imag
    data = np.array([omegas.real*HARTREE2EV, p_ave])
    np.savetxt('test_0157_bse_h2o_rhf_cis_pyscf.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0157_bse_h2o_rhf_cis_pyscf.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))
    
    nao_td  = bse_iter(mf=gto_mf, gto=mol, verbosity=0, xc_code='CIS')

    polariz = -nao_td.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*HARTREE2EV, polariz])
    np.savetxt('test_0157_bse_h2o_rhf_cis_nao.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0157_bse_h2o_rhf_cis_nao.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=5e-5, rtol=5e-2), 
      msg="{}".format(abs(data_ref-data).sum()/data.size ) ) 

def main(args):
    # dataset
    testset = get_datasets(args)
    batch_size = len(testset)

    amp = args.deg * math.pi / 180.0
    w = torch.randn(batch_size, 3)
    w = w / w.norm(p=2, dim=1, keepdim=True) * amp
    t = torch.rand(batch_size, 3) * args.max_trans

    if args.format == 'wv':
        # the output: twist vectors.
        R = ptlk.so3.exp(w) # (N, 3) --> (N, 3, 3)
        G = torch.zeros(batch_size, 4, 4)
        G[:, 3, 3] = 1
        G[:, 0:3, 0:3] = R
        G[:, 0:3, 3] = t

        x = ptlk.se3.log(G) # --> (N, 6)
    else:
        # rotation-vector and translation-vector
        x = torch.cat((w, t), dim=1)

    numpy.savetxt(args.outfile, x, delimiter=',') 

def test_0163_bse_h2o_spin2_uhf_rpa(self):
    """ This  """
    mol = gto.M(verbose=1,atom='O 0 0 0; H 0 0.489 1.074; H 0 0.489 -1.074',basis='cc-pvdz',spin=2)

    gto_mf = scf.UKS(mol)
    gto_mf.xc = 'hf'
    gto_mf.kernel()
    gto_td = tddft.dRPA(gto_mf)
    gto_td.nstates = 190
    gto_td.kernel()

    omegas = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    p_ave = -polariz_freq_osc_strength(gto_td.e, gto_td.oscillator_strength(), omegas).imag
    data = np.array([omegas.real*HARTREE2EV, p_ave])
    np.savetxt('test_0163_bse_h2o_spin2_uhf_rpa_pyscf.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0163_bse_h2o_spin2_uhf_rpa_pyscf.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))
    
    nao_td  = bse_iter(mf=gto_mf, gto=mol, verbosity=0, xc_code='RPA')

    polariz = -nao_td.comp_polariz_inter_ave(omegas).imag
    data = np.array([omegas.real*HARTREE2EV, polariz])
    np.savetxt('test_0163_bse_h2o_spin2_uhf_rpa_nao.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0163_bse_h2o_spin2_uhf_rpa_nao.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-5, rtol=1e-3), \
      msg="{}".format(abs(data_ref-data).sum()/data.size)) 

def test_tddft_upkc(self):
    """ This is a comparison of two equivalent ways of computing the polarizability for water molecule """
    td = tddft_iter_2ord(label='water', cd=dirname(abspath(__file__)),jcutoff=7,xc_code='RPA',verbosity=0)
    omegas = np.arange(0.0,1.0,0.005)+1j*0.01

    pxx1 = -td.comp_polariz_nonin_xx(omegas).imag
    data1 = np.array([omegas.real*27.2114, pxx1])
    np.savetxt('water.tddft_iter_nonin.txt', data1.T)
    #print('    td.rf0_ncalls ', td.rf0_ncalls)
    #print(' td.matvec_ncalls ', td.matvec_ncalls)


    pxx1 = -td.comp_polariz_inter_xx(omegas).imag
    data1 = np.array([omegas.real*27.2114, pxx1])
    np.savetxt('water.tddft_iter_unit.txt', data1.T)
    #print('    td.rf0_ncalls ', td.rf0_ncalls)
    #print(' td.matvec_ncalls ', td.matvec_ncalls)

    pxx2 = td.polariz_upkc(omegas)
    wp = np.zeros((2*pxx2.shape[1]+1,pxx2.shape[0]))
    wp[0,:] = omegas.real*27.2114
    wp[1:pxx2.shape[1]+1,:] = pxx2.real.T
    wp[pxx2.shape[1]+1:,:] = pxx2.imag.T
    np.savetxt('water.tddft_iter_upkc.txt', wp.T)
    #print('    td.rf0_ncalls ', td.rf0_ncalls)
    #print(' td.matvec_ncalls ', td.matvec_ncalls)

    pxx3 = -td.polariz_dckcd(omegas).imag
    data1 = np.array([omegas.real*27.2114, pxx3])
    np.savetxt('water.tddft_iter_dckcd.txt', data1.T)
    #print('    td.rf0_ncalls ', td.rf0_ncalls)
    #print(' td.matvec_ncalls ', td.matvec_ncalls) 

def test_138_h2o_rks_rpa_pb(self):
    """ This  """
    mol = gto.M(verbose=1,atom='O 0 0 0; H 0 0.489 1.074; H 0 0.489 -1.074',basis='cc-pvdz')
    gto_mf = dft.RKS(mol)
    gto_mf.kernel()
    nao_mf = tddft_iter(gto=mol, mf=gto_mf)
    comega = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    pnonin = -nao_mf.comp_polariz_inter_ave(comega, verbosity=0).imag
    data = np.array([comega.real*HARTREE2EV, pnonin])
    np.savetxt('test_138_h2o_rks_pz_pb.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_138_h2o_rks_pz_pb.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3)) 

def test_156_bse_h2b_uks_nonin(self):
    """ This  example is with discrepancies... """
    mol = gto.M(verbose=1,atom='B 0 0 0; H 0 0.489 1.074; H 0 0.489 -1.074',basis='cc-pvdz',spin=3)

    gto_mf = scf.UHF(mol)
    gto_mf.kernel()

    omegas = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    p_ave = -polariz_nonin_ave(gto_mf, mol, omegas).imag
    data = np.array([omegas.real*HARTREE2EV, p_ave])
    np.savetxt('test_0156_bse_h2b_uhf_nonin_pyscf.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0156_bse_h2b_uhf_nonin_pyscf.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))
    
    nao_td  = bse_iter(mf=gto_mf, gto=mol, verbosity=0, xc_code='RPA')

    polariz = -nao_td.polariz_nonin_ave_matelem(omegas).imag
    data = np.array([omegas.real*HARTREE2EV, polariz])
    np.savetxt('test_0156_bse_h2b_uhf_nonin_matelem.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0156_bse_h2b_uhf_nonin_matelem.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3))
    
    polariz = -nao_td.comp_polariz_nonin_ave(omegas).imag
    data = np.array([omegas.real*HARTREE2EV, polariz])
    np.savetxt('test_0156_bse_h2b_uhf_nonin_nao.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_0156_bse_h2b_uhf_nonin_nao.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3)) 

def test_132_h2o_uhf_nonin_pb(self):
    """ This  """
    mol = gto.M(verbose=1,atom='O 0 0 0; H 0 0.489 1.074; H 0 0.489 -1.074',basis='cc-pvdz',spin=0)
    gto_mf = scf.UHF(mol)
    gto_mf.kernel()
    nao_mf = chi0_matvec(gto=mol, mf=gto_mf, tol_loc=1e-5, tol_biloc=1e-7)
    comega = np.arange(0.0, 2.0, 0.01) + 1j*0.03
    pnonin = -nao_mf.comp_polariz_nonin_ave(comega, verbosity=1).imag
    data = np.array([comega.real*HARTREE2EV, pnonin])
    np.savetxt('test_132_h2o_uhf_nonin_pb.txt', data.T, fmt=['%f','%f'])
    data_ref = np.loadtxt('test_132_h2o_uhf_nonin_pb.txt-ref').T
    self.assertTrue(np.allclose(data_ref, data, atol=1e-6, rtol=1e-3)) 

def test_bse_gto_vs_nao_nonin_0081(self):
    """ Non-interacting case """
    #print(__name__, 'gto.mo_energy', gto_mf.mo_energy)
    nao_td = bse_iter(mf=gto_mf, gto=mol, verbosity=0, perform_gw=True)
    
    omegas = np.linspace(0.0,2.0,450)+1j*0.04
    p_iter = -nao_td.comp_polariz_nonin_ave(omegas).imag
    data = np.array([omegas.real*27.2114, p_iter])
    np.savetxt('be.bse_iter.omega.nonin.ave.txt', data.T, fmt=['%f','%f']) 

def test_0066_al_atom(self):
    """ Spin-resolved case GW procedure. """
    gw = gw_c(mf=gto_mf, gto=mol, verbosity=0, niter_max_ev=16, nocc=3, nvrt=3)
    self.assertEqual(gw.nspin, 2)
    gw.kernel_gw()
    #gw.report()
    np.savetxt('eigvals_g0w0_pyscf_rescf_al_0066.txt', gw.mo_energy_gw[0,:,:].T)
    #ev_ref = np.loadtxt('eigvals_g0w0_pyscf_rescf_al_0066.txt-ref').T
    #for n2e,n2r in zip(gw.mo_energy_gw[0], ev_ref):
    #  for e,r in zip(n2e,n2r): self.assertAlmostEqual(e, r)