Construct a ``CavityAnalytic`` instance from the Pynac lattice element
497410
def from_dynacRepr(cls, pynacRepr):
        
        cavID = int(pynacRepr[1][0][0])
        xesln = float(pynacRepr[1][1][0])
        phase = float(pynacRepr[1][1][1])
        fieldReduction = float(pynacRepr[1][1][2])
        isec = int(pynacRepr[1][1][3])
        return cls(phase, fieldReduction, cavID, xesln, isec)
725003
def as_dict(self):
        
        d = {"vasp_version": self.vasp_version,
             "has_vasp_completed": self.converged,
             "nsites": len(self.final_structure)}
        comp = self.final_structure.composition
        d["unit_cell_formula"] = comp.as_dict()
        d["reduced_cell_formula"] = Composition(comp.reduced_formula).as_dict()
        d["pretty_formula"] = comp.reduced_formula
        symbols = [s.split()[1] for s in self.potcar_symbols]
        symbols = [re.split(r"_", s)[0] for s in symbols]
        d["is_hubbard"] = self.is_hubbard
        d["hubbards"] = self.hubbards

        unique_symbols = sorted(list(set(self.atomic_symbols)))
        d["elements"] = unique_symbols
        d["nelements"] = len(unique_symbols)

        d["run_type"] = self.run_type

        vin = {"incar": {k: v for k, v in self.incar.items()},
               "crystal": self.initial_structure.as_dict(),
               "kpoints": self.kpoints.as_dict()}
        actual_kpts = [{"abc": list(self.actual_kpoints[i]),
                        "weight": self.actual_kpoints_weights[i]}
                       for i in range(len(self.actual_kpoints))]
        vin["kpoints"]["actual_points"] = actual_kpts
        vin["nkpoints"] = len(actual_kpts)
        vin["potcar"] = [s.split(" ")[1] for s in self.potcar_symbols]
        vin["potcar_spec"] = self.potcar_spec
        vin["potcar_type"] = [s.split(" ")[0] for s in self.potcar_symbols]
        vin["parameters"] = {k: v for k, v in self.parameters.items()}
        vin["lattice_rec"] = self.lattice_rec.as_dict()
        d["input"] = vin

        nsites = len(self.final_structure)

        try:
            vout = {"ionic_steps": self.ionic_steps,
                    "final_energy": self.final_energy,
                    "final_energy_per_atom": self.final_energy / nsites,
                    "crystal": self.final_structure.as_dict(),
                    "efermi": self.efermi}
        except (ArithmeticError, TypeError):
            vout = {"ionic_steps": self.ionic_steps,
                    "final_energy": self.final_energy,
                    "final_energy_per_atom": None,
                    "crystal": self.final_structure.as_dict(),
                    "efermi": self.efermi}

        if self.eigenvalues:
            eigen = {str(spin): v.tolist()
                     for spin, v in self.eigenvalues.items()}
            vout["eigenvalues"] = eigen
            (gap, cbm, vbm, is_direct) = self.eigenvalue_band_properties
            vout.update(dict(bandgap=gap, cbm=cbm, vbm=vbm,
                             is_gap_direct=is_direct))

            if self.projected_eigenvalues:
                vout['projected_eigenvalues'] = {
                    str(spin): v.tolist()
                    for spin, v in self.projected_eigenvalues.items()}

        vout['epsilon_static'] = self.epsilon_static
        vout['epsilon_static_wolfe'] = self.epsilon_static_wolfe
        vout['epsilon_ionic'] = self.epsilon_ionic
        d['output'] = vout
        return jsanitize(d, strict=True)
829192
def _initialise_approximator(self, closest_point=None,
        wavelengths_required=None, rescale=True, **kwargs):
        

        if self._initialised and not kwargs.get("force", False):
            logger.debug("Ignoring call to re-initialise approximator because "
                "we already have.")
            return self._subset_bounds


        logger.info("Initialising approximator near {0}".format(closest_point))

        N = self.grid_points.size
        dtype = [(name, '<f8') for name in self.grid_points.dtype.names]
        grid_points = self.grid_points.astype(dtype).view(float).reshape(N, -1)

        # If closest_point is given then we will slice a small part of the grid.
        if closest_point is not None:

            default = 1.0
            grid_subset = self._configuration.get("settings", 
                { "grid_subset": default }).get("grid_subset", default)

            # Use nearest N points in interpolator.
            # If grid subset is a fraction, scale it to real numbers.
            if 1 >= grid_subset > 0:
                grid_subset = int(np.round(grid_subset * N))

            logger.debug("Using {0} nearest points for interpolator".format(
                grid_subset))

            # Use closest N points.
            distance = np.sum(np.abs(grid_points - closest_point)/
                np.ptp(grid_points, axis=0), axis=1)
            grid_indices = np.argsort(distance)[:grid_subset]

        else:
            grid_indices = np.ones(N, dtype=bool)

        # What wavelengths will be required? Assume all if not specified.
        if wavelengths_required is not None:
            mask = np.zeros(self.wavelengths.size, dtype=bool)
            for start, end in wavelengths_required:
                idx = np.clip(
                    self.wavelengths.searchsorted([start, end]) + [0, 1],
                    0, self.wavelengths.size)
                mask[idx[0]:idx[1]] = True
        else:
            mask = np.ones(self.wavelengths.size, dtype=bool)

        # Apply model masks.
        mask *= self._model_mask()

        # Slice small part of the intensities grid.
        intensities = np.memmap(
            self._configuration["model_grid"]["intensities"],
            dtype="float32", mode="r", shape=(N, self.wavelengths.size))
        subset = np.copy(intensities[grid_indices, :])
        subset[:, ~mask] = np.nan
        del intensities

        # Create an interpolator.
        try:
            interpolator = interpolate.LinearNDInterpolator(
                grid_points[grid_indices], subset, rescale=rescale)

        except TypeError:
            logger.warn("Could not rescale the LinearNDInterpolator because "\
                "you need a newer version of scipy")
            interpolator = interpolate.LinearNDInterpolator(
                grid_points[grid_indices], subset)

        generate.init()
        generate.wavelengths.append(self.wavelengths)
        generate.intensities.append(interpolator)
        
        self._initialised = True

        # Return the subset boundaries of the grid.
        self._subset_bounds = {}
        for name in self.grid_points.dtype.names:
            points = self.grid_points[name][grid_indices]
            self._subset_bounds[name] = (min(points), max(points))

        return self._subset_bounds
446621
def _construct_tensorflow_expression(slvr_cfg, feed_data, device, shard):
    
    zero = tf.constant(0)
    src_count = zero
    src_ph_vars = feed_data.src_ph_vars

    LSA = feed_data.local

    polarisation_type = slvr_cfg['polarisation_type']

    # Pull RIME inputs out of the feed staging_area
    # of the relevant shard, adding the feed once
    # inputs to the dictionary
    D = LSA.feed_many[shard].get_to_attrdict()
    D.update({k: fo.var for k, fo in LSA.feed_once.iteritems()})

    with tf.device(device):
        # Infer chunk dimensions
        model_vis_shape = tf.shape(D.model_vis)
        ntime, nbl, nchan, npol = [model_vis_shape[i] for i in range(4)]

        # Infer float and complex type
        FT, CT = D.uvw.dtype, D.model_vis.dtype

        # Compute sine and cosine of parallactic angles
        # for the beam
        beam_sin, beam_cos = rime.parallactic_angle_sin_cos(
                                        D.parallactic_angles)

        # Compute sine and cosine of feed rotation angle
        feed_sin, feed_cos = rime.parallactic_angle_sin_cos(
                                        D.parallactic_angles[:, :] +
                                        D.feed_angles[None, :])

        # Compute feed rotation
        feed_rotation = rime.feed_rotation(feed_sin, feed_cos, CT=CT,
                                           feed_type=polarisation_type)

    def antenna_jones(lm, stokes, alpha, ref_freq):
        """
        Compute the jones terms for each antenna.

        lm, stokes and alpha are the source variables.
        """

        # Compute the complex phase
        cplx_phase = rime.phase(lm, D.uvw, D.frequency, CT=CT)

        # Check for nans/infs in the complex phase
        phase_msg = ("Check that '1 - l**2  - m**2 >= 0' holds "
                    "for all your lm coordinates. This is required "
                    "for 'n = sqrt(1 - l**2 - m**2) - 1' "
                    "to be finite.")

        phase_real = tf.check_numerics(tf.real(cplx_phase), phase_msg)
        phase_imag = tf.check_numerics(tf.imag(cplx_phase), phase_msg)

        # Compute the square root of the brightness matrix
        # (as well as the sign)
        bsqrt, sgn_brightness = rime.b_sqrt(stokes, alpha,
            D.frequency, ref_freq, CT=CT,
            polarisation_type=polarisation_type)

        # Check for nans/infs in the bsqrt
        bsqrt_msg = ("Check that your stokes parameters "
                    "satisfy I**2 >= Q**2 + U**2 + V**2. "
                    "Montblanc performs a cholesky decomposition "
                    "of the brightness matrix and the above must "
                    "hold for this to produce valid values.")

        bsqrt_real = tf.check_numerics(tf.real(bsqrt), bsqrt_msg)
        bsqrt_imag = tf.check_numerics(tf.imag(bsqrt), bsqrt_msg)

        # Compute the direction dependent effects from the beam
        ejones = rime.e_beam(lm, D.frequency,
            D.pointing_errors, D.antenna_scaling,
            beam_sin, beam_cos,
            D.beam_extents, D.beam_freq_map, D.ebeam)

        deps = [phase_real, phase_imag, bsqrt_real, bsqrt_imag]
        deps = [] # Do nothing for now

        # Combine the brightness square root, complex phase,
        # feed rotation and beam dde's
        with tf.control_dependencies(deps):
            antenna_jones = rime.create_antenna_jones(bsqrt, cplx_phase,
                                                    feed_rotation, ejones, FT=FT)
            return antenna_jones, sgn_brightness

    # While loop condition for each point source type
    def point_cond(coherencies, npsrc, src_count):
        return tf.less(npsrc, src_ph_vars.npsrc)

    def gaussian_cond(coherencies, ngsrc, src_count):
        return tf.less(ngsrc, src_ph_vars.ngsrc)

    def sersic_cond(coherencies, nssrc, src_count):
        return tf.less(nssrc, src_ph_vars.nssrc)

    # While loop bodies
    def point_body(coherencies, npsrc, src_count):
        """ Accumulate visiblities for point source batch """
        S = LSA.sources['npsrc'][shard].get_to_attrdict()

        # Maintain source counts
        nsrc = tf.shape(S.point_lm)[0]
        src_count += nsrc
        npsrc +=  nsrc

        ant_jones, sgn_brightness = antenna_jones(S.point_lm,
            S.point_stokes, S.point_alpha, S.point_ref_freq)
        shape = tf.ones(shape=[nsrc,ntime,nbl,nchan], dtype=FT)
        coherencies = rime.sum_coherencies(D.antenna1, D.antenna2,
            shape, ant_jones, sgn_brightness, coherencies)

        return coherencies, npsrc, src_count

    def gaussian_body(coherencies, ngsrc, src_count):
        """ Accumulate coherencies for gaussian source batch """
        S = LSA.sources['ngsrc'][shard].get_to_attrdict()

        # Maintain source counts
        nsrc = tf.shape(S.gaussian_lm)[0]
        src_count += nsrc
        ngsrc += nsrc

        ant_jones, sgn_brightness = antenna_jones(S.gaussian_lm,
            S.gaussian_stokes, S.gaussian_alpha, S.gaussian_ref_freq)
        gauss_shape = rime.gauss_shape(D.uvw, D.antenna1, D.antenna2,
            D.frequency, S.gaussian_shape)
        coherencies = rime.sum_coherencies(D.antenna1, D.antenna2,
            gauss_shape, ant_jones, sgn_brightness, coherencies)

        return coherencies, ngsrc, src_count

    def sersic_body(coherencies, nssrc, src_count):
        """ Accumulate coherencies for sersic source batch """
        S = LSA.sources['nssrc'][shard].get_to_attrdict()

        # Maintain source counts
        nsrc = tf.shape(S.sersic_lm)[0]
        src_count += nsrc
        nssrc += nsrc

        ant_jones, sgn_brightness = antenna_jones(S.sersic_lm,
            S.sersic_stokes, S.sersic_alpha, S.sersic_ref_freq)
        sersic_shape = rime.sersic_shape(D.uvw, D.antenna1, D.antenna2,
            D.frequency, S.sersic_shape)
        coherencies = rime.sum_coherencies(D.antenna1, D.antenna2,
            sersic_shape, ant_jones, sgn_brightness, coherencies)

        return coherencies, nssrc, src_count

    with tf.device(device):
        base_coherencies = tf.zeros(shape=[ntime,nbl,nchan,npol], dtype=CT)

        # Evaluate point sources
        summed_coherencies, npsrc, src_count = tf.while_loop(
            point_cond, point_body,
            [base_coherencies, zero, src_count])

        # Evaluate gaussians
        summed_coherencies, ngsrc, src_count = tf.while_loop(
            gaussian_cond, gaussian_body,
            [summed_coherencies, zero, src_count])

        # Evaluate sersics
        summed_coherencies, nssrc, src_count = tf.while_loop(
            sersic_cond, sersic_body,
            [summed_coherencies, zero, src_count])

        # Post process visibilities to produce model visibilites and chi squared
        model_vis, chi_squared = rime.post_process_visibilities(
            D.antenna1, D.antenna2, D.direction_independent_effects, D.flag,
            D.weight, D.model_vis, summed_coherencies, D.observed_vis)

    # Create enstaging_area operation
    put_op = LSA.output.put_from_list([D.descriptor, model_vis, chi_squared])

    # Return descriptor and enstaging_area operation
    return D.descriptor, put_op
683946
def build_boolCoeff(self):
        ''' 
        '''
        # coefficients for hill functions from boolean update rules
        self.boolCoeff = collections.OrderedDict([(s,[]) for s in self.varNames.keys()])
        # parents
        self.pas = collections.OrderedDict([(s,[]) for s in self.varNames.keys()])
        #
        for key in self.boolRules.keys():
            rule = self.boolRules[key]
            self.pas[key] = self.parents_from_boolRule(rule)
            pasIndices = [self.varNames[pa] for pa in self.pas[key]]
            # check whether there are coupling matrix entries for each parent
            for g in range(self.dim):
                if g in pasIndices:
                    if np.abs(self.Coupl[self.varNames[key],g]) < 1e-10:
                        raise ValueError('specify coupling value for '+str(key)+' <- '+str(g))
                else:
                    if np.abs(self.Coupl[self.varNames[key],g]) > 1e-10:
                        raise ValueError('there should be no coupling value for '+str(key)+' <- '+str(g))
            if self.verbosity > 1:
                settings.m(0,'...'+key)
                settings.m(0,rule)
                settings.m(0,rule_pa)
            # now evaluate coefficients
            for tuple in list(itertools.product([False,True],repeat=len(self.pas[key]))):
                if self.process_rule(rule,self.pas[key],tuple):
                    self.boolCoeff[key].append(tuple)
            #
            if self.verbosity > 1:
                settings.m(0,self.boolCoeff[key])
1270774
def perform_es_corr(self, lattice, q, step=1e-4):
        
        logger.info("Running Freysoldt 2011 PC calculation (should be " "equivalent to sxdefectalign)")
        logger.debug("defect lattice constants are (in angstroms)" + str(lattice.abc))

        [a1, a2, a3] = ang_to_bohr * np.array(lattice.get_cartesian_coords(1))
        logging.debug("In atomic units, lat consts are (in bohr):" + str([a1, a2, a3]))
        vol = np.dot(a1, np.cross(a2, a3))  # vol in bohr^3

        def e_iso(encut):
            gcut = eV_to_k(encut)  # gcut is in units of 1/A
            return scipy.integrate.quad(lambda g: self.q_model.rho_rec(g * g)**2, step, gcut)[0] * (q**2) / np.pi

        def e_per(encut):
            eper = 0
            for g2 in generate_reciprocal_vectors_squared(a1, a2, a3, encut):
                eper += (self.q_model.rho_rec(g2)**2) / g2
            eper *= (q**2) * 2 * round(np.pi, 6) / vol
            eper += (q**2) * 4 * round(np.pi, 6) \
                * self.q_model.rho_rec_limit0 / vol
            return eper

        eiso = converge(e_iso, 5, self.madetol, self.energy_cutoff)
        logger.debug("Eisolated : %f", round(eiso, 5))

        eper = converge(e_per, 5, self.madetol, self.energy_cutoff)

        logger.info("Eperiodic : %f hartree", round(eper, 5))
        logger.info("difference (periodic-iso) is %f hartree", round(eper - eiso, 6))
        logger.info("difference in (eV) is %f", round((eper - eiso) * hart_to_ev, 4))

        es_corr = round((eiso - eper) / self.dielectric * hart_to_ev, 6)
        logger.info("Defect Correction without alignment %f (eV): ", es_corr)
        return es_corr
761524
def get_x_objs(self, *dynac_type):
        
        return [i for i in self.lattice for y in dynac_type if dynac_from_ele(i) == y]
199037
def _get_lattice_parameters(lattice):
    

    return np.array(np.sqrt(np.dot(lattice.T, lattice).diagonal()),
                    dtype='double')
1259648
def _build(self):
    
    if "w" not in self._initializers:
      stddev = 1 / math.sqrt(np.prod(self._shape))
      self._initializers["w"] = tf.truncated_normal_initializer(stddev=stddev)

    self._w = tf.get_variable("w",
                              shape=self._shape,
                              dtype=self._dtype,
                              initializer=self._initializers["w"],
                              partitioner=self._partitioners.get("w", None),
                              regularizer=self._regularizers.get("w", None))
    return self._w
1069007
def _build_model(self):
        

        # Set up LSTM modules
        self.lstms = nn.ModuleList(
            [
                RNN(
                    num_classes=0,
                    num_tokens=self.word_dict.s,
                    emb_size=self.settings["emb_dim"],
                    lstm_hidden=self.settings["hidden_dim"],
                    attention=self.settings["attention"],
                    dropout=self.settings["dropout"],
                    bidirectional=self.settings["bidirectional"],
                    use_cuda=self.settings["host_device"] in self._gpu,
                )
            ]
            * self.settings["relation_arity"]
        )

        if "input_dim" not in self.settings:
            raise ValueError("Model parameter input_dim cannot be None.")

        # Set up final linear layer
        self.sparse_linear = SparseLinear(
            self.settings["input_dim"], self.cardinality, self.settings["bias"]
        )
775106
def make_fc(dim_in, hidden_dim, use_gn=False):
    '''
        
    '''
    if use_gn:
        fc = nn.Linear(dim_in, hidden_dim, bias=False)
        nn.init.kaiming_uniform_(fc.weight, a=1)
        return nn.Sequential(fc, group_norm(hidden_dim))
    fc = nn.Linear(dim_in, hidden_dim)
    nn.init.kaiming_uniform_(fc.weight, a=1)
    nn.init.constant_(fc.bias, 0)
    return fc
341185
def cublasZgemm(handle, transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc):
    

    status = _libcublas.cublasZgemm_v2(handle,
                                       _CUBLAS_OP[transa],
                                       _CUBLAS_OP[transb], m, n, k, 
                                       ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                         alpha.imag)),
                                       int(A), lda, int(B), ldb,
                                       ctypes.byref(cuda.cuDoubleComplex(beta.real,
                                                                         beta.imag)),
                                       int(C), ldc)
    cublasCheckStatus(status)
494489
def generator(ngf, nc, no_bias=True, fix_gamma=True, eps=1e-5 + 1e-12, z_dim=100, activation='sigmoid'):
    '''
    '''
    BatchNorm = mx.sym.BatchNorm
    rand = mx.sym.Variable('rand')

    rand = mx.sym.Reshape(rand, shape=(-1, z_dim, 1, 1))

    g1 = mx.sym.Deconvolution(rand, name='gen1', kernel=(5,5), stride=(2,2),target_shape=(2,2), num_filter=ngf*8, no_bias=no_bias)
    gbn1 = BatchNorm(g1, name='genbn1', fix_gamma=fix_gamma, eps=eps)
    gact1 = mx.sym.Activation(gbn1, name="genact1", act_type="relu")

    g2 = mx.sym.Deconvolution(gact1, name='gen2', kernel=(5,5), stride=(2,2),target_shape=(4,4), num_filter=ngf*4, no_bias=no_bias)
    gbn2 = BatchNorm(g2, name='genbn2', fix_gamma=fix_gamma, eps=eps)
    gact2 = mx.sym.Activation(gbn2, name='genact2', act_type='relu')

    g3 = mx.sym.Deconvolution(gact2, name='gen3', kernel=(5,5), stride=(2,2), target_shape=(8,8), num_filter=ngf*2, no_bias=no_bias)
    gbn3 = BatchNorm(g3, name='genbn3', fix_gamma=fix_gamma, eps=eps)
    gact3 = mx.sym.Activation(gbn3, name='genact3', act_type='relu')

    g4 = mx.sym.Deconvolution(gact3, name='gen4', kernel=(5,5), stride=(2,2), target_shape=(16,16), num_filter=ngf, no_bias=no_bias)
    gbn4 = BatchNorm(g4, name='genbn4', fix_gamma=fix_gamma, eps=eps)
    gact4 = mx.sym.Activation(gbn4, name='genact4', act_type='relu')

    g5 = mx.sym.Deconvolution(gact4, name='gen5', kernel=(5,5), stride=(2,2), target_shape=(32,32), num_filter=nc, no_bias=no_bias)
    gout = mx.sym.Activation(g5, name='genact5', act_type=activation)

    return gout
803175
def Construct(self):  # pylint: disable-msg=C0103
        """"""
        # Parse the GDML
        self.world = self.gdml_parser.GetWorldVolume()

        # Create sensitive detector
        self.sensitive_detector = ScintSD()

        # Get logical volume for X view, then attach SD
        my_lv = G4.G4LogicalVolumeStore.GetInstance().GetVolumeID(1)
        assert my_lv.GetName() == "ScintillatorBarX"
        my_lv.SetSensitiveDetector(self.sensitive_detector)

        # Get logical volume for Y view, then attach SD
        my_lv = G4.G4LogicalVolumeStore.GetInstance().GetVolumeID(2)
        assert my_lv.GetName() == "ScintillatorBarY"
        my_lv.SetSensitiveDetector(self.sensitive_detector)

        my_lv = G4.G4LogicalVolumeStore.GetInstance().GetVolumeID(0)
        assert my_lv.GetName() == "SteelPlane"

        # field
        self.field_manager = G4.G4FieldManager()
        self.my_field = MagneticField.WandsToroidField(self.field_polarity)
        self.field_manager.SetDetectorField(self.my_field)
        self.field_manager.CreateChordFinder(self.my_field)
        my_lv.SetFieldManager(self.field_manager, False)

        self.log.info("Materials:")
        self.log.info(G4.G4Material.GetMaterialTable())

        # Return pointer to world volume
        return self.world
392031
def create(input_width, input_height, input_channels=1, output_dim=512):
    
    def instantiate(**_):
        return NatureCnn(
            input_width=input_width, input_height=input_height, input_channels=input_channels,
            output_dim=output_dim
        )

    return ModelFactory.generic(instantiate)
1325181
private void calcInterfacesCrystal(StructureInterfaceList set, double cutoff) {


		// initialising debugging vars
		long start = -1;
		long end = -1;
		int trialCount = 0;
		int skippedRedundant = 0;
		int skippedAUsNoOverlap = 0;
		int skippedChainsNoOverlap = 0;
		int skippedSelfEquivalent = 0;

		// The bounding boxes of all AUs of the unit cell
		UnitCellBoundingBox bbGrid = new UnitCellBoundingBox(numOperatorsSg, numPolyChainsAu);;
		// we calculate all the bounds of each of the asym units, those will then be reused and translated
		bbGrid.setBbs(structure, ops, INCLUDE_HETATOMS);


		// if not crystallographic there's no search to do in other cells, only chains within "AU" will be checked
		if (!searchBeyondAU) numCells = 0;

		boolean verbose = logger.isDebugEnabled();

		if (verbose) {
			trialCount = 0;
			start= System.currentTimeMillis();
			int neighbors = (2*numCells+1)*(2*numCells+1)*(2*numCells+1)-1;
			int auTrials = (numPolyChainsAu*(numPolyChainsAu-1))/2;
			int trials = numPolyChainsAu*numOperatorsSg*numPolyChainsAu*neighbors;
			logger.debug("Chain clash trials within original AU: "+auTrials);
			logger.debug(
					"Chain clash trials between the original AU and the neighbouring "+neighbors+
					" whole unit cells ("+numCells+" neighbours)" +
					"(2x"+numPolyChainsAu+"chains x "+numOperatorsSg+"AUs x "+neighbors+"cells) : "+trials);
			logger.debug("Total trials: "+(auTrials+trials));
		}

		List<Chain> polyChains = structure.getPolyChains();

		for (int a=-numCells;a<=numCells;a++) {
			for (int b=-numCells;b<=numCells;b++) {
				for (int c=-numCells;c<=numCells;c++) {

					Point3i trans = new Point3i(a,b,c);
					Vector3d transOrth = new Vector3d(a,b,c);
					if (a!=0 || b!=0 || c!=0) {
						// we avoid doing the transformation for 0,0,0 (in case it's not crystallographic)
						this.crystallographicInfo.getCrystalCell().transfToOrthonormal(transOrth);
					}

					UnitCellBoundingBox bbGridTrans = bbGrid.getTranslatedBbs(transOrth);

					for (int n=0;n<numOperatorsSg;n++) {

						// short-cut strategies
						// 1) we skip first of all if the bounding boxes of the AUs don't overlap
						if (!bbGrid.getAuBoundingBox(0).overlaps(bbGridTrans.getAuBoundingBox(n), cutoff)) {
							skippedAUsNoOverlap++;
							continue;
						}

						// 2) we check if we didn't already see its equivalent symmetry operator partner
						CrystalTransform tt = new CrystalTransform(this.crystallographicInfo.getSpaceGroup(), n);
						tt.translate(trans);
						if (isRedundantTransform(tt)) {
							skippedRedundant++;
							continue;
						}
						addVisitedTransform(tt);


						boolean selfEquivalent = false;

						// 3) an operator can be "self redundant" if it is the inverse of itself (involutory, e.g. all pure 2-folds with no translation)
						if (tt.isEquivalent(tt)) {
							logger.debug("Transform {} is equivalent to itself, will skip half of i-chains to j-chains comparisons", tt.toString());
							// in this case we can't skip the operator, but we can skip half of the matrix comparisons e.g. j>i
							// we set a flag and do that within the loop below
							selfEquivalent = true;
						}

						StringBuilder builder = null;
						if (verbose) builder = new StringBuilder(String.valueOf(tt)).append(" ");

						// Now that we know that boxes overlap and operator is not redundant, we have to go to the details
						int contactsFound = 0;

						for (int j=0;j<numPolyChainsAu;j++) {

							for (int i=0;i<numPolyChainsAu;i++) { // we only have to compare the original asymmetric unit to every full cell around

								if(selfEquivalent && (j>i)) {
									// in case of self equivalency of the operator we can safely skip half of the matrix
									skippedSelfEquivalent++;
									continue;
								}
								// special case of original AU, we don't compare a chain to itself
								if (n==0 && a==0 && b==0 && c==0 && i==j) continue;

								// before calculating the AtomContactSet we check for overlap, then we save putting atoms into the grid
								if (!bbGrid.getChainBoundingBox(0,i).overlaps(bbGridTrans.getChainBoundingBox(n,j),cutoff)) {
									skippedChainsNoOverlap++;
									if (verbose) {
										builder.append(".");
									}
									continue;
								}

								trialCount++;

								// finally we've gone through all short-cuts and the 2 chains seem to be close enough:
								// we do the calculation of contacts
								Chain chaini = polyChains.get(i);
								Chain chainj = polyChains.get(j);

								if (n!=0 || a!=0 || b!=0 || c!=0) {
									Matrix4d mJCryst = new Matrix4d(ops[n]);
									translate(mJCryst, transOrth);
									chainj = (Chain)chainj.clone();
									Calc.transform(chainj,mJCryst);
								}

								StructureInterface interf = calcContacts(chaini, chainj, cutoff, tt, builder);
								if (interf == null) {
									continue;
								}

								contactsFound++;
								if(this.hasNcsOps()) {
									StructureInterface interfNcsRef = findNcsRef(interf);
									set.addNcsEquivalent(interf,interfNcsRef);
								} else {
									set.add(interf);
								}
							}
						}

						if( verbose ) {
							if (a==0 && b==0 && c==0 && n==0)
								builder.append(" "+contactsFound+"("+(numPolyChainsAu*(numPolyChainsAu-1))/2+")");
							else if (selfEquivalent)
								builder.append(" "+contactsFound+"("+(numPolyChainsAu*(numPolyChainsAu+1))/2+")");
							else
								builder.append(" "+contactsFound+"("+numPolyChainsAu*numPolyChainsAu+")");

							logger.debug(builder.toString());
						}
					}
				}
			}
		}

		end = System.currentTimeMillis();
		logger.debug("\n"+trialCount+" chain-chain clash trials done. Time "+(end-start)/1000+"s");
		logger.debug("  skipped (not overlapping AUs)       : "+skippedAUsNoOverlap);
		logger.debug("  skipped (not overlapping chains)    : "+skippedChainsNoOverlap);
		logger.debug("  skipped (sym redundant op pairs)    : "+skippedRedundant);
		logger.debug("  skipped (sym redundant self op)     : "+skippedSelfEquivalent);
		logger.debug("Found "+set.size()+" interfaces.");
	}
991078
def from_dynacRepr(cls, pynacRepr):
        
        L = float(pynacRepr[1][0][0])
        B = float(pynacRepr[1][0][1])
        aperRadius = float(pynacRepr[1][0][2])
        return cls(L, B, aperRadius)