BSE calculations with Haydock can be restarted only if we have the
        excitonic Hamiltonian and the HAYDR_SAVE file.
517979
def restart(self):
        
        # TODO: This version seems to work but the main output file is truncated
        # TODO: Handle restart if CG method is used
        # TODO: restart should receive a list of critical events
        # the log file is complete though.
        irdvars = {}

        # Move the BSE blocks to indata.
        # This is done only once at the end of the first run.
        # Successive restarts will use the BSR|BSC files in the indir directory
        # to initialize the excitonic Hamiltonian
        count = 0
        for ext in ("BSR", "BSC"):
            ofile = self.outdir.has_abiext(ext)
            if ofile:
                count += 1
                irdvars.update(irdvars_for_ext(ext))
                self.out_to_in(ofile)

        if not count:
            # outdir does not contain the BSR|BSC file.
            # This means that num_restart > 1 and the files should be in task.indir
            count = 0
            for ext in ("BSR", "BSC"):
                ifile = self.indir.has_abiext(ext)
                if ifile:
                    count += 1

            if not count:
                raise self.RestartError("%s: Cannot find BSR|BSC files in %s" % (self, self.indir))

        # Rename HAYDR_SAVE files
        count = 0
        for ext in ("HAYDR_SAVE", "HAYDC_SAVE"):
            ofile = self.outdir.has_abiext(ext)
            if ofile:
                count += 1
                irdvars.update(irdvars_for_ext(ext))
                self.out_to_in(ofile)

        if not count:
            raise self.RestartError("%s: Cannot find the HAYDR_SAVE file to restart from." % self)

        # Add the appropriate variable for restarting.
        self.set_vars(irdvars)

        # Now we can resubmit the job.
        #self.history.info("Will restart from %s", restart_file)
        return self._restart()
676261
def inter_spin_hamiltonian(self, u_int, J_coup):
        """"""
        J_coup *= u_int
        h_int  = (u_int - 2*J_coup)/2.*self.oper['sumSz2']
        h_int += J_coup*self.oper['sumSz-sp2']
        h_int -= J_coup/2.*self.oper['sumSz-or2']
        h_int -= J_coup*self.oper['Sfliphop']

        return h_int
1217436
def _iso_ic_on_monochromatic(
        electron_energy, seed_energy, seed_edensity, gamma_energy
    ):
        
        photE0 = (seed_energy / mec2).decompose().value
        phn = seed_edensity

        # electron_energy = electron_energy[:, None]
        gamma_energy = gamma_energy[:, None]
        photE0 = photE0[:, None, None]
        phn = phn[:, None, None]

        b = 4 * photE0 * electron_energy
        w = gamma_energy / electron_energy
        q = w / (b * (1 - w))
        fic = (
            2 * q * np.log(q)
            + (1 + 2 * q) * (1 - q)
            + (1.0 / 2.0) * (b * q) ** 2 * (1 - q) / (1 + b * q)
        )

        gamint = (
            fic
            * heaviside(1 - q)
            * heaviside(q - 1.0 / (4 * electron_energy ** 2))
        )
        gamint[np.isnan(gamint)] = 0.0

        if phn.size > 1:
            phn = phn.to(1 / (mec2_unit * u.cm ** 3)).value
            gamint = trapz_loglog(gamint * phn / photE0, photE0, axis=0)  # 1/s
        else:
            phn = phn.to(mec2_unit / u.cm ** 3).value
            gamint *= phn / photE0 ** 2
            gamint = gamint.squeeze()

        # gamint /= mec2.to('erg').value

        # r0 = (e**2 / m_e / c**2).to('cm')
        # sigt = ((8 * np.pi) / 3 * r0**2).cgs
        sigt = 6.652458734983284e-25
        c = 29979245800.0

        gamint *= (3.0 / 4.0) * sigt * c / electron_energy ** 2

        return gamint
524664
def reset(self):
    
    self._poolingActivation = numpy.zeros((self._numColumns), dtype="int32")
    self._poolingColumns = []
    self._overlapDutyCycles = numpy.zeros(self._numColumns, dtype=realDType)
    self._activeDutyCycles = numpy.zeros(self._numColumns, dtype=realDType)
    self._minOverlapDutyCycles = numpy.zeros(self._numColumns,
                                             dtype=realDType)
    self._minActiveDutyCycles = numpy.zeros(self._numColumns,
                                            dtype=realDType)
    self._boostFactors = numpy.ones(self._numColumns, dtype=realDType)
776553
def cublasDsyrk(handle, uplo, trans, n, k, alpha, A, lda, beta, C, ldc):
    
    
    status = _libcublas.cublasDsyrk_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       n, k, ctypes.byref(cuda.cuFloatComplex(alpha.real,     
                                                                              alpha.imag)),
                                       int(A), lda, 
                                       ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                        beta.imag)),
                                       int(C), ldc)
    cublasCheckStatus(status)
789269
def reset(self):
        
        self.solved = False
        self.niter = 0
        self.iter_mis = []
        self.F = None
        self.system.dae.factorize = True
200715
private void rotateShiftAtoms(Atom[] ca){

		for (int i  = 0 ; i < ca.length; i++){
			Atom c = ca[i];

			Calc.rotate(c,currentRotMatrix);
			Calc.shift(c,currentTranMatrix);
			//System.out.println("after " + c);
			ca[i] = c;
		}
		//System.out.println("after " + ca[0]);
	}
1339582
def save(self, calc, session):
        '''
        
        '''
        checksum = calc.get_checksum()

        try:
            existing_calc = session.query(model.Calculation).filter(model.Calculation.checksum == checksum).one()
        except NoResultFound:
            pass
        else:
            del calc
            return None, "This calculation already exists!"

        if not calc.download_size:
            for f in calc.related_files:
                calc.download_size += os.stat(f).st_size

        ormcalc = model.Calculation(checksum = checksum)

        if calc._calcset:
            ormcalc.meta_data = model.Metadata(chemical_formula = calc.info['standard'], download_size = calc.download_size)

            for child in session.query(model.Calculation).filter(model.Calculation.checksum.in_(calc._calcset)).all():
                ormcalc.children.append(child)
            ormcalc.siblings_count = len(ormcalc.children)
            ormcalc.nested_depth = calc._nested_depth

        else:
            # prepare phonon data for saving
            # this is actually a dict to list conversion TODO re-structure this
            if calc.phonons['modes']:
                phonons_json = []

                for bzpoint, frqset in calc.phonons['modes'].items():
                    # re-orientate eigenvectors
                    for i in range(0, len(calc.phonons['ph_eigvecs'][bzpoint])):
                        for j in range(0, len(calc.phonons['ph_eigvecs'][bzpoint][i])//3):
                            eigv = array([calc.phonons['ph_eigvecs'][bzpoint][i][j*3], calc.phonons['ph_eigvecs'][bzpoint][i][j*3+1], calc.phonons['ph_eigvecs'][bzpoint][i][j*3+2]])
                            R = dot( eigv, calc.structures[-1].cell ).tolist()
                            calc.phonons['ph_eigvecs'][bzpoint][i][j*3], calc.phonons['ph_eigvecs'][bzpoint][i][j*3+1], calc.phonons['ph_eigvecs'][bzpoint][i][j*3+2] = [round(x, 3) for x in R]

                    try: irreps = calc.phonons['irreps'][bzpoint]
                    except KeyError:
                        empty = []
                        for i in range(len(frqset)):
                            empty.append('')
                        irreps = empty

                    phonons_json.append({  'bzpoint':bzpoint, 'freqs':frqset, 'irreps':irreps, 'ph_eigvecs':calc.phonons['ph_eigvecs'][bzpoint]  })
                    if bzpoint == '0 0 0':
                        phonons_json[-1]['ir_active'] = calc.phonons['ir_active']
                        phonons_json[-1]['raman_active'] = calc.phonons['raman_active']
                    if calc.phonons['ph_k_degeneracy']:
                        phonons_json[-1]['ph_k_degeneracy'] = calc.phonons['ph_k_degeneracy'][bzpoint]

                ormcalc.phonons = model.Phonons()
                ormcalc.spectra.append( model.Spectra(kind = model.Spectra.PHONON, eigenvalues = json.dumps(phonons_json)) )

            # prepare electron data for saving TODO re-structure this
            for task in ['dos', 'bands']: # projected?
                if calc.electrons[task]:
                    calc.electrons[task] = calc.electrons[task].todict()

            if calc.electrons['dos'] or calc.electrons['bands']:
                ormcalc.electrons = model.Electrons(gap = calc.info['bandgap'])
                if 'bandgaptype' in calc.info:
                    ormcalc.electrons.is_direct = 1 if calc.info['bandgaptype'] == 'direct' else -1
                ormcalc.spectra.append(model.Spectra(
                    kind = model.Spectra.ELECTRON,
                    dos = json.dumps(calc.electrons['dos']),
                    bands = json.dumps(calc.electrons['bands']),
                    projected = json.dumps(calc.electrons['projected']),
                    eigenvalues = json.dumps(calc.electrons['eigvals'])
                ))

            # construct ORM for other props
            calc.related_files = list(map(virtualize_path, calc.related_files))
            ormcalc.meta_data = model.Metadata(location = calc.info['location'], finished = calc.info['finished'], raw_input = calc.info['input'], modeling_time = calc.info['duration'], chemical_formula = html_formula(calc.info['standard']), download_size = calc.download_size, filenames = json.dumps(calc.related_files))

            codefamily = model.Codefamily.as_unique(session, content = calc.info['framework'])
            codeversion = model.Codeversion.as_unique(session, content = calc.info['prog'])

            codeversion.instances.append( ormcalc.meta_data )
            codefamily.versions.append( codeversion )

            pot = model.Pottype.as_unique(session, name = calc.info['H'])
            pot.instances.append(ormcalc)
            ormcalc.recipinteg = model.Recipinteg(kgrid = calc.info['k'], kshift = calc.info['kshift'], smearing = calc.info['smear'], smeartype = calc.info['smeartype'])
            ormcalc.basis = model.Basis(kind = calc.info['ansatz'], content = json.dumps(calc.electrons['basis_set']) if calc.electrons['basis_set'] else None)
            ormcalc.energy = model.Energy(convergence = json.dumps(calc.convergence), total = calc.info['energy'])

            ormcalc.spacegroup = model.Spacegroup(n=calc.info['ng'])
            ormcalc.struct_ratios = model.Struct_ratios(chemical_formula=calc.info['standard'], formula_units=calc.info['expanded'], nelem=calc.info['nelem'], dimensions=calc.info['dims'])
            if len(calc.tresholds) > 1:
                ormcalc.struct_optimisation = model.Struct_optimisation(tresholds=json.dumps(calc.tresholds), ncycles=json.dumps(calc.ncycles))

            for n, ase_repr in enumerate(calc.structures):
                is_final = True if n == len(calc.structures)-1 else False
                struct = model.Structure(step = n, final = is_final)

                s = cell_to_cellpar(ase_repr.cell)
                struct.lattice = model.Lattice(a=s[0], b=s[1], c=s[2], alpha=s[3], beta=s[4], gamma=s[5], a11=ase_repr.cell[0][0], a12=ase_repr.cell[0][1], a13=ase_repr.cell[0][2], a21=ase_repr.cell[1][0], a22=ase_repr.cell[1][1], a23=ase_repr.cell[1][2], a31=ase_repr.cell[2][0], a32=ase_repr.cell[2][1], a33=ase_repr.cell[2][2])

                #rmts =      ase_repr.get_array('rmts') if 'rmts' in ase_repr.arrays else [None for j in range(len(ase_repr))]
                charges =   ase_repr.get_array('charges') if 'charges' in ase_repr.arrays else [None for j in range(len(ase_repr))]
                magmoms =   ase_repr.get_array('magmoms') if 'magmoms' in ase_repr.arrays else [None for j in range(len(ase_repr))]
                for n, i in enumerate(ase_repr):
                    struct.atoms.append( model.Atom( number=chemical_symbols.index(i.symbol), x=i.x, y=i.y, z=i.z, charge=charges[n], magmom=magmoms[n] ) )

                ormcalc.structures.append(struct)
            # TODO Forces

        ormcalc.uigrid = model.Grid(info=json.dumps(calc.info))

        # tags ORM
        uitopics = []
        for entity in self.hierarchy:

            if not entity['creates_topic']:
                continue

            if entity['multiple'] or calc._calcset:
                for item in calc.info.get( entity['source'], [] ):
                    uitopics.append( model.topic(cid=entity['cid'], topic=item) )
            else:
                topic = calc.info.get(entity['source'])
                if topic or not entity['optional']:
                    uitopics.append( model.topic(cid=entity['cid'], topic=topic) )

        uitopics = [model.Topic.as_unique(session, cid=x.cid, topic="%s" % x.topic) for x in uitopics]

        ormcalc.uitopics.extend(uitopics)

        if calc._calcset:
            session.add(ormcalc)
        else:
            session.add_all([codefamily, codeversion, pot, ormcalc])

        session.commit()
        del calc, ormcalc
        return checksum, None
1108129
private boolean isSp3Carbon(IAtom atom, int deg) {
        Integer elem = atom.getAtomicNumber();
        Integer hcnt = atom.getImplicitHydrogenCount();
        if (elem == null || hcnt == null) return false;
        if (elem == 6 && hcnt <= 1 && deg + hcnt == 4) {
            // more expensive check, look one out and see if we have any
            // duplicate terminal neighbors
            List<IAtom> terminals = new ArrayList<>();
            for (IBond bond : container.getConnectedBondsList(atom)) {
                IAtom nbr = bond.getOther(atom);
                if (container.getConnectedBondsCount(nbr) == 1) {
                    for (IAtom terminal : terminals) {
                        if (Objects.equals(terminal.getAtomicNumber(),
                                           nbr.getAtomicNumber()) &&
                            Objects.equals(terminal.getMassNumber(),
                                           nbr.getMassNumber()) &&
                            Objects.equals(terminal.getFormalCharge(),
                                           nbr.getFormalCharge()) &&
                            Objects.equals(terminal.getImplicitHydrogenCount(),
                                           nbr.getImplicitHydrogenCount())) {
                            return false;
                        }
                    }
                    terminals.add(nbr);
                }
            }
            return true;
        }
        return false;
    }
819731
def solve_unchecked(self, sense=None):
        
        if sense is not None:
            self.set_objective_sense(sense)

        parm = swiglpk.glp_smcp()
        swiglpk.glp_init_smcp(parm)
        if self._do_presolve:
            parm.presolve = swiglpk.GLP_ON
            self._do_presolve = False
        else:
            parm.presolve = swiglpk.GLP_OFF

        parm.tol_bnd = self._feasibility_tolerance
        parm.tol_dj = self._optimality_tolerance

        logger.debug('Solving using glp_simplex()')
        r = swiglpk.glp_simplex(self._p, parm)
        if r in (swiglpk.GLP_ENOPFS, swiglpk.GLP_ENODFS):
            self._result = Result(self, r)
            return self._result
        elif r != 0:
            raise GLPKError('glp_simplex: {}'.format(r))

        if swiglpk.glp_get_num_int(self._p) == 0:
            self._result = Result(self)
        else:
            # The intopt MILP solver needs an optimal solution to the LP
            # relaxation. Therefore, glp_simplex has to run before glp_intopt
            # for MILP problems.
            logger.debug('Solving using glp_intopt()')
            parm = swiglpk.glp_iocp()
            swiglpk.glp_init_iocp(parm)
            parm.tol_int = self._integrality_tolerance
            r = swiglpk.glp_intopt(self._p, parm)
            if r != 0:
                raise GLPKError('glp_intopt: {}'.format(r))

            self._result = MIPResult(self)

        return self._result
1355412
def _cleanup(self, ref):
        ''
        self.transitor_states[self._weakref_holder[ref]] -= 1
        del self._weakref_holder[ref]
17926
public static void rebuildAFPChain(AFPChain afpChain, Atom[] ca1, Atom[] ca2){

		if ( afpChain.getAlgorithmName() == null) {
			afpChain.setAlgorithmName(DEFAULT_ALGORITHM_NAME);
		}
		if ( afpChain.getVersion() == null){
			afpChain.setVersion("1.0");
		}

		int blockNum  = afpChain.getBlockNum();
		int ca1Length = afpChain.getCa1Length();
		int ca2Length = afpChain.getCa2Length();

		int minLength = Math.min(ca1Length, ca2Length);
		int[][][] optAln = new int[blockNum][2][minLength];

		int[][][] blockResList = afpChain.getBlockResList();
		if ( blockResList == null){
			blockResList = new int[blockNum][2][minLength];
		}
		int[] optLen = afpChain.getOptLen();

		String[][][] pdbAln = afpChain.getPdbAln();
		int[] verifiedOptLen = null;
		if ( optLen != null)
		  verifiedOptLen = afpChain.getOptLen().clone();
		else {
			logger.warn("did not find optimal alignment, building up empty alignment.");
			optLen = new int[1];
			optLen[0] = 0;
		}
		for (int blockNr = 0 ; blockNr < blockNum ; blockNr++){

			//System.out.println("got block " + blockNr + " size: " + optLen[blockNr]);
			int verifiedEQR = -1;
			for ( int eqrNr = 0 ; eqrNr < optLen[blockNr] ; eqrNr++ ){
				String pdbResnum1 = pdbAln[blockNr][0][eqrNr];
				String pdbResnum2 = pdbAln[blockNr][1][eqrNr];

				//System.out.println(blockNr + " " + eqrNr + " got resnum: " + pdbResnum1 + " " + pdbResnum2);
				String[] spl1 = pdbResnum1.split(":");
				String[] spl2 = pdbResnum2.split(":");

				String chain1 = spl1[0];
				String pdbres1 = spl1[1];

				String chain2 = spl2[0];
				String pdbres2 = spl2[1];

				int pos1 = getPositionForPDBresunm(pdbres1,chain1,ca1);
				int pos2 = getPositionForPDBresunm(pdbres2,chain2,ca2);

				if ( pos1 == -1 || pos2 == -1 ){
					// this can happen when parsing old files that contained Calcium atoms...
					logger.warn("pos1: {} (residue {}), pos2: {} (residue {}), should never be -1. Probably parsing an old file.",
							pos1, pdbResnum1, pos2, pdbResnum2);
					verifiedOptLen[blockNr]-- ;
					continue;
				}

				verifiedEQR++;
				//System.out.println(blockNr + " " + eqrNr + " " + pos1 + " " + pos2);
				optAln[blockNr][0][verifiedEQR] = pos1;
				optAln[blockNr][1][verifiedEQR] = pos2;
				blockResList[blockNr][0][verifiedEQR] = pos1;
				blockResList[blockNr][1][verifiedEQR] = pos2;
			}
		}

		afpChain.setOptLen(verifiedOptLen);
		afpChain.setOptAln(optAln);
		afpChain.setBlockResList(blockResList);
		// build up alignment image:
		AFPAlignmentDisplay.getAlign(afpChain, ca1, ca2);


	}
305018
def QA_SU_save_stock_terminated(client=DATABASE):
    '''
    
    '''

    # 🛠todo 已经失效从wind 资讯里获取
    # 这个函数已经失效
    print("！！！ tushare 这个函数已经失效！！！")
    df = QATs.get_terminated()
    #df = QATs.get_suspended()
    print(
        " Get stock terminated from tushare,stock count is %d  (终止上市股票列表)" %
        len(df)
    )
    coll = client.stock_terminated
    client.drop_collection(coll)
    json_data = json.loads(df.reset_index().to_json(orient='records'))
    coll.insert(json_data)
    print(" 保存终止上市股票列表 到 stock_terminated collection， OK")
699384
double getBetheFreeEnergy() {
        // 
        // G_{Bethe} = \sum_a \sum_{x_a} - b(x_a) ln \chi(x_a)
        //              + \sum_a \sum_{x_a} b(x_a) ln b(x_a)
        //              + \sum_i (n_i - 1) \sum_{x_i} b(x_i) ln b(x_i)
        //           = \sum_a \sum_{x_a} b(x_a) ln (b(x_a) / \chi(x_a))
        //              + \sum_i (n_i - 1) \sum_{x_i} b(x_i) ln b(x_i)
        //
        //     where n_i is the number of neighbors of the variable x_i,
        //     b(x_a) and b(x_i) are normalized distributions and x_a is 
        //     the set of variables participating in factor a. 
        //
                
        double bethe = 0.0;
        for (int a=0; a<fg.getFactors().size(); a++) {
            Factor f = fg.getFactors().get(a);
            if (!(f instanceof GlobalFactor)) {
                int numConfigs = f.getVars().calcNumConfigs();
                VarTensor beliefs = getFactorBeliefs(a);
                for (int c=0; c<numConfigs; c++) {                
                    // Since we want multiplication by 0 to always give 0 (not the case for Double.POSITIVE_INFINITY or Double.NaN.
                    double b_c = beliefs.getValue(c);
                    if (b_c != s.zero()) {
                        double r_b_c = s.toReal(b_c);
                        double log_b_c = s.toLogProb(b_c);
                        double log_chi_c = f.getLogUnormalizedScore(c);
                        bethe += r_b_c * (log_b_c - log_chi_c);
                    }
                }
            } else {
                VarTensor[] inMsgs = getMsgs(f.getId(), msgs, IN_MSG);
                bethe += ((GlobalFactor) f).getExpectedLogBelief(inMsgs);
            }
        }
        for (int v=0; v<fg.getVars().size(); v++) {
            Var var = fg.getVars().get(v);
            int numNeighbors = bg.numNbsT1(v);
            VarTensor beliefs = getVarBeliefs(v);
            double sum = 0.0;
            for (int c=0; c<var.getNumStates(); c++) {
                double b_c = beliefs.getValue(c);
                if (b_c != s.zero()) {
                    double r_b_c = s.toReal(b_c);
                    double log_b_c = s.toLogProb(b_c);
                    sum += r_b_c * log_b_c;
                }
            }
            bethe -= (numNeighbors - 1) * sum;
        }
        
        assert !Double.isNaN(bethe);        
        return bethe;
    }
160567
def eta_from_seebeck(seeb,Lambda):
    
    from scipy.optimize import fsolve
    out = fsolve(lambda x: (seebeck_spb(x,Lambda) - abs(seeb)) ** 2, 1.,full_output=True)
    return out[0][0]
51246
def _fix_lsm_bitspersample(self, parent):
        
        if self.code != 258 or self.count != 2:
            return
        # TODO: test this case; need example file
        log.warning('TiffTag %i: correcting LSM bitspersample tag', self.code)
        value = struct.pack('<HH', *self.value)
        self.valueoffset = struct.unpack('<I', value)[0]
        parent.filehandle.seek(self.valueoffset)
        self.value = struct.unpack('<HH', parent.filehandle.read(4))
576819
def selfconsistency(self, u_int, J_coup, mean_field_prev=None):
        """"""
        if mean_field_prev is None:
            mean_field_prev = np.array([self.param['ekin']]*2)
        hlog = [mean_field_prev]

        self.oper['Hint'] = self.inter_spin_hamiltonian(u_int, J_coup)
        converging = True
        half_fill = (self.param['populations'] == 0.5).all()
        while converging:
            if half_fill:
                self.update_H(hlog[-1], self.param['lambda'])
            else:
                res = root(self.restriction, self.param['lambda'], (hlog[-1]))#, method='lm')
                if not res.success:
                    res.x = res.x * 0.5 + 0.5*self.param['lambda']
                    self.update_H(self.mean_field()*0.5 + 0.5*hlog[-1], res.x)
                    print('fail', self.param['populations'][3:5])
                if (self.quasiparticle_weight() < 0.001).all():
                    return hlog
                self.param['lambda'] = res.x

            hlog.append(self.mean_field())
            converging = (abs(hlog[-1] - hlog[-2]) > self.param['tol']).all() \
                    or (abs(self.restriction(self.param['lambda'], hlog[-1])) > self.param['tol']).all()

        return hlog