def get_trajectory_datasets(tgrp, *fields):
    '''Return a list of new/existing datasets corresponding to the given fields

       **Arguments:**

       tgrp
            The trajectory group

       fields
            A list of fields, i.e. pairs of name and row_shape.
    '''
    result = []
    for name, row_shape in fields:
        if name in tgrp:
            ds = tgrp[name]
            if ds.shape[1:] != row_shape:
                raise TypeError('The shape of the existing dataset is not compatible with the new data.')
            if log.do_medium:
                log('Found an existing dataset %s in group %s with %i rows.' % (name, tgrp.name, ds.shape[0]))
        else:
            if log.do_high:
                log('Creating new dataset %s with row shape %s' % (name, row_shape))
            # Create a new dataset
            shape = (0,) + row_shape
            maxshape = (None,) + row_shape
            ds = tgrp.create_dataset(name, shape, maxshape=maxshape, dtype=float)
        result.append(ds)
    return result

def g09log_to_hdf5(f, fn_log):
    """Convert Gaussian09 BOMD log file to Yaff HDF5 format.

       **Arguments:**

       f
            An open and writable HDF5 file.

       fn_log
            The name of the Gaussian log file.
    """
    with log.section('G09H5'):
        if log.do_medium:
            log('Loading Gaussian 09 file \'%s\' into \'trajectory\' of HDF5 file \'%s\'' % (
                fn_log, f.filename
            ))

        # First make sure the HDF5 file has a system description that is consistent
        # with the XYZ file.
        if 'system' not in f:
            raise ValueError('The HDF5 file must contain a system group.')
        if 'numbers' not in f['system']:
            raise ValueError('The HDF5 file must have a system group with atomic numbers.')
        natom = f['system/numbers'].shape[0]

        # Take care of the trajectory group
        tgrp = get_trajectory_group(f)

        # Take care of the pos and vel datasets
        dss = get_trajectory_datasets(tgrp,
            ('pos', (natom, 3)),
            ('vel', (natom, 3)),
            ('frc', (natom, 3)),
            ('time', (1,)),
            ('step', (1,)),
            ('epot', (1,)),
            ('ekin', (1,)),
            ('etot', (1,)),
        )
        ds_pos, ds_vel, ds_frc, ds_time, ds_step, ds_epot, ds_ekin, ds_etot = dss

        # Load frame by frame
        row = get_last_trajectory_row(dss)
        for numbers, pos, vel, frc, time, step, epot, ekin, etot in _iter_frames_g09(fn_log):
            if (numbers != f['system/numbers']).any():
                log.warn('The element numbers of the HDF5 and LOG file do not match.')
            write_to_dataset(ds_pos, pos, row)
            write_to_dataset(ds_vel, vel, row)
            write_to_dataset(ds_frc, frc, row)
            write_to_dataset(ds_time, time, row)
            write_to_dataset(ds_step, step, row)
            write_to_dataset(ds_epot, epot, row)
            write_to_dataset(ds_ekin, ekin, row)
            write_to_dataset(ds_etot, etot, row)
            row += 1

        # Check number of rows
        check_trajectory_rows(tgrp, dss, row)

 def __call__(self, iterative):
     if log.do_medium:
         if self.time0 is None:
             self.time0 = time.time()
             if log.do_medium:
                 log.hline()
                 log('counter  Walltime')
                 log.hline()
         log('%7i %10.1f' % (
             iterative.counter,
             time.time() - self.time0,
         ))

    def apply(self, par_table, cpar_table, scale_table, mixing_rules, system, ff_args):
        # Prepare the atomic parameters
        amps = np.zeros(system.nffatype, float)
        bs = np.zeros(system.nffatype, float)
        for i in xrange(system.nffatype):
            key = (system.ffatypes[i],)
            par_list = par_table.get(key, [])
            if len(par_list) == 0:
                if log.do_warning:
                    log.warn('No EXPREP parameters found for ffatype %s.' % system.ffatypes[i])
            else:
                amps[i], bs[i] = par_list[0]

        # Prepare the cross parameters
        amp_cross = np.zeros((system.nffatype, system.nffatype), float)
        b_cross = np.zeros((system.nffatype, system.nffatype), float)
        for i0 in xrange(system.nffatype):
            for i1 in xrange(i0+1):
                cpar_list = cpar_table.get((system.ffatypes[i0], system.ffatypes[i1]), [])
                if len(cpar_list) == 0:
                    if log.do_high:
                        log('No EXPREP cross parameters found for ffatypes %s,%s. Mixing rule will be used' % (system.ffatypes[i0], system.ffatypes[i1]))
                else:
                    amp_cross[i0,i1], b_cross[i0,i1] = cpar_list[0]
                    if i0 != i1:
                        amp_cross[i1,i0], b_cross[i1,i0] = cpar_list[0]

        # Prepare the global parameters
        scalings = Scalings(system, scale_table[1], scale_table[2], scale_table[3])
        amp_mix, amp_mix_coeff = mixing_rules['A']
        if amp_mix == 0:
            amp_mix_coeff = 0.0
        elif amp_mix == 1:
            amp_mix_coeff = amp_mix_coeff[0]
        b_mix, b_mix_coeff = mixing_rules['B']
        if b_mix == 0:
            b_mix_coeff = 0.0
        elif b_mix == 1:
            b_mix_coeff = b_mix_coeff[0]

        # Get the part. It should not exist yet.
        part_pair = ff_args.get_part_pair(PairPotExpRep)
        if part_pair is not None:
            raise RuntimeError('Internal inconsistency: the EXPREP part should not be present yet.')

        pair_pot = PairPotExpRep(
            system.ffatype_ids, amp_cross, b_cross, ff_args.rcut, ff_args.tr,
            amps, amp_mix, amp_mix_coeff, bs, b_mix, b_mix_coeff,
        )
        nlist = ff_args.get_nlist(system)
        part_pair = ForcePartPair(system, nlist, scalings, pair_pot)
        ff_args.parts.append(part_pair)

    def update(self, dx, dg):
        tmp = dg - np.dot(self.hessian, dx)

        denom = np.dot(tmp, dx)
        if abs(denom) > 1e-5 * np.linalg.norm(dx) * np.linalg.norm(tmp):
            if log.do_debug:
                log("Updating SR1 Hessian.       denom=%10.3e" % denom)
            self.hessian += np.outer(tmp, tmp) / denom
            return True
        else:
            if log.do_high:
                log("Skipping SR1 update because denom=%10.3e is not big enough." % denom)
            return False

    def make_step(self):
        # get relevant hessian information
        evals, evecs = self.hessian.get_spectrum()
        if log.do_high:
            log(" lowest eigen value: %7.1e" % evals.min())
            log("highest eigen value: %7.1e" % evals.max())
        # convert gradient to eigen basis
        grad_eigen = np.dot(evecs.T, self.g_old)

        while True:
            # Find the step with the given radius. If the hessian is positive
            # definite and the unconstrained step is smaller than the trust
            # radius, this step is returned
            delta_eigen = solve_trust_radius(grad_eigen, evals, self.trust_radius)
            radius = np.linalg.norm(delta_eigen)

            # convert the step to user basis
            delta_x = np.dot(evecs, delta_eigen)

            # compute the function and gradient at the new position
            x = self.x_old + delta_x
            f, g = self.fun(x, True)

            # compute the change in function value
            delta_f = f - self.f_old
            # compute the change in norm of the gradient
            delta_norm_g = np.linalg.norm(g) - np.linalg.norm(self.g_old)
            # must_shrink is a parameter to control the trust radius
            must_shrink = False

            if delta_f > 0:
                # The function must decrease, if not the trust radius is too big.
                if log.do_high:
                    log("Function increases.")
                must_shrink = True

            if self.trust_radius < self.small_radius and delta_norm_g > 0:
                # When the trust radius becomes small, the numerical noise on
                # the energy may be too large to detect an increase energy.
                # In that case the norm of the gradient is used instead.
                if log.do_high:
                    log("Gradient norm increases.")
                must_shrink = True

            if must_shrink:
                self.trust_radius *= 0.5
                while self.trust_radius >= radius:
                    self.trust_radius *= 0.5
                if self.trust_radius < self.too_small_radius:
                    raise RuntimeError("The trust radius becomes too small. Is the potential energy surface smooth?")
            else:
                # If we get here, we are done with the trust radius loop.
                if log.do_high:
                    log.hline()
                # It is fine to increase the trust radius a little after a
                # successful step.
                if self.trust_radius < self.initial_trust_radius:
                    self.trust_radius *= 2.0
                # Return the results of the successful step
                return x, f, g

def get_trajectory_group(f):
    '''Create or return an existing trajectory group

       **Arguments:**

       f
            An open HDF5 File or Group object.
    '''
    if 'trajectory' not in f:
        if log.do_high:
            log('Creating new trajectory datagroup in %s.' % f.filename)
        tgrp = f.create_group('trajectory')
    else:
        tgrp = f['trajectory']
        if log.do_high:
            log('Using existing trajectory datagroup in %s.' % f.filename)
    return tgrp

    def _verify_hooks(self):
        with log.section('ENSEM'):
            thermo = None
            index_thermo = 0
            baro = None
            index_baro = 0

            # Look for the presence of a thermostat and/or barostat
            if hasattr(self.hooks, '__len__'):
                for index, hook in enumerate(self.hooks):
                    if hook.method == 'thermostat':
                        thermo = hook
                        index_thermo = index
                    elif hook.method == 'barostat':
                        baro = hook
                        index_baro = index
            elif self.hooks is not None:
                if self.hooks.method == 'thermostat':
                    thermo = self.hooks
                elif self.hooks.method == 'barostat':
                    baro = self.hooks

            # If both are present, delete them and generate TBCombination element
            if thermo is not None and baro is not None:
                from yaff.sampling.npt import TBCombination
                if log.do_warning:
                    log.warn('Both thermostat and barostat are present separately and will be merged')
                del self.hooks[max(index_thermo, index_thermo)]
                del self.hooks[min(index_thermo, index_baro)]
                self.hooks.append(TBCombination(thermo, baro))

            if hasattr(self.hooks, '__len__'):
                for hook in self.hooks:
                    if hook.name == 'TBCombination':
                        thermo = hook.thermostat
                        baro = hook.barostat
            elif self.hooks is not None:
                if self.hooks.name == 'TBCombination':
                    thermo = self.hooks.thermostat
                    baro = self.hooks.barostat

            if log.do_warning:
                if thermo is not None:
                    log('Temperature coupling achieved through ' + str(thermo.name) + ' thermostat')
                if baro is not None:
                    log('Pressure coupling achieved through ' + str(baro.name) + ' barostat')

 def propagate(self):
     # Update the Hessian
     assert not self.g is self.g_old
     assert not self.x is self.x_old
     hessian_safe = self.hessian.update(self.x - self.x_old, self.g - self.g_old)
     if not hessian_safe:
         # Reset the Hessian completely
         if log.do_high:
             log("Resetting hessian due to failed update.")
         self.hessian = SR1HessianModel(len(self.x))
         self.trust_radius = self.initial_trust_radius
     # Move new to old
     self.x_old = self.x
     self.f_old = self.f
     self.g_old = self.g
     # Compute a step
     self.x, self.f, self.g = self.make_step()
     return BaseOptimizer.propagate(self)

 def __call__(self, iterative):
     if log.do_medium:
         if self.time0 is None:
             self.time0 = time.time()
             if log.do_medium:
                 log.hline()
                 log('Cons.Err. =&the root of the ratio of the variance on the conserved quantity and the variance on the kinetic energy.')
                 log('d-rmsd    =&the root-mean-square displacement of the atoms.')
                 log('g-rmsd    =&the root-mean-square gradient of the energy.')
                 log('counter  Cons.Err.       Temp     d-RMSD     g-RMSD   Walltime')
                 log.hline()
         log('%7i %10.5f %s %s %s %10.1f' % (
             iterative.counter,
             iterative.cons_err,
             log.temperature(iterative.temp),
             log.length(iterative.rmsd_delta),
             log.force(iterative.rmsd_gpos),
             time.time() - self.time0,
         ))

 def __call__(self, iterative):
     if log.do_medium:
         if self.time0 is None:
             self.time0 = time.time()
             if log.do_medium:
                 log.hline()
                 log('Conv.val. =&the highest ratio of a convergence criterion over its threshold.')
                 log('N         =&the number of convergence criteria that is not met.')
                 log('Worst     =&the name of the convergence criterion that is worst.')
                 log('counter  Conv.val.  N           Worst     Energy   Walltime')
                 log.hline()
         log('%7i % 10.3e %2i %15s %s %10.1f' % (
             iterative.counter,
             iterative.dof.conv_val,
             iterative.dof.conv_count,
             iterative.dof.conv_worst,
             log.energy(iterative.epot),
             time.time() - self.time0,
         ))

    def from_hdf5(cls, f):
        '''Create a system from an HDF5 file/group containing a system group

           **Arguments:**

           f
                An open h5.File object with a system group. The system group
                must at least contain a numbers and pos dataset.
        '''
        sgrp = f['system']
        kwargs = {
            'numbers': sgrp['numbers'][:],
            'pos': sgrp['pos'][:],
        }
        for key in 'scopes', 'scope_ids', 'ffatypes', 'ffatype_ids', 'bonds', 'rvecs', 'charges', 'masses':
            if key in sgrp:
                kwargs[key] = sgrp[key][:]
        if log.do_high:
            log('Read system parameters from %s.' % f.filename)
        return cls(**kwargs)

    def __init__(self, system, comlist=None):
        '''
           Parameters
           ----------

           system
                An instance of the ``System`` class.
           comlist
                An optional layer to derive centers of mass from the atomic positions.
                These centers of mass are used as input for the first layer, the relative
                vectors.
        '''
        ForcePart.__init__(self, 'valence', system)
        self.comlist = comlist
        self.dlist = DeltaList(system if comlist is None else comlist)
        self.iclist = InternalCoordinateList(self.dlist)
        self.vlist = ValenceList(self.iclist)
        if log.do_medium:
            with log.section('FPINIT'):
                log('Force part: %s' % self.name)
                log.hline()

 def _init_log(self):
     if log.do_medium:
         log('Unit cell')
         log.hline()
         log('Number of periodic dimensions: %i' % self.cell.nvec)
         lengths, angles = self.cell.parameters
         names = 'abc'
         for i in xrange(len(lengths)):
             log('Cell parameter %5s: %10s' % (names[i], log.length(lengths[i])))
         names = 'alpha', 'beta', 'gamma'
         for i in xrange(len(angles)):
             log('Cell parameter %5s: %10s' % (names[i], log.angle(angles[i])))
         log.hline()
         log.blank()

 def _init_derived_scopes(self):
     if self.scope_ids is None:
         if len(self.scopes) != self.natom:
             raise TypeError('When the scope_ids are derived automatically, the length of the scopes list must match the number of atoms.')
         lookup = {}
         scopes = []
         self.scope_ids = np.zeros(self.natom, int)
         for i in xrange(self.natom):
             scope = self.scopes[i]
             scope_id = lookup.get(scope)
             if scope_id is None:
                 scope_id = len(scopes)
                 scopes.append(scope)
                 lookup[scope] = scope_id
             self.scope_ids[i] = scope_id
         self.scopes = scopes
     for scope in self.scopes:
         check_name(scope)
     # check the range of the ids
     if self.scope_ids.min() != 0 or self.scope_ids.max() != len(self.scopes)-1:
         raise ValueError('The ffatype_ids have incorrect bounds.')
     if log.do_medium:
         log('The following scopes are present in the system:')
         log.hline()
         log('                 Scope   ID   Number of atoms')
         log.hline()
         for scope_id, scope in enumerate(self.scopes):
             log('%22s  %3i       %3i' % (scope, scope_id, (self.scope_ids==scope_id).sum()))
         log.hline()
         log.blank()