Chimie Théorique » QMX

import sys, os

class Qmx(Calculator):

        self._constraints=None

    def get_energy_subsystem(self, path, calculator, atoms, force_consistent):

        # go to directory and calculate energies

        # perform energy calculations

        e_sys_lo = self.get_energy_subsystem("system.low-level", self.calculator_low_system, embed.get_system(), force_consistent)

        e_cl_lo  = self.get_energy_subsystem("cluster.low-level", self.calculator_low_cluster, embed.get_cluster(), force_consistent)

        # calculate energies

        energy = e_sys_lo - e_cl_lo + e_cl_hi

        # print energies

        # set energies and return

        if force_consistent:

            self.energy_free = energy

        # get forces for the three systems

        f_sys_lo = self.get_forces_subsystem("system.low-level", self.calculator_low_system, embed.get_system())

        f_cl_lo  = self.get_forces_subsystem("cluster.low-level", self.calculator_low_cluster, embed.get_cluster())

        # forces correction for the atoms

        f_cl = f_cl_hi - f_cl_lo

        # lo-sys + (hi-lo)

            iat_cl = atom_map_sys_cl[iat_sys]

            if iat_cl > -1:

                f_sys_lo[iat_sys] += f_cl[iat_cl]

        # correct gradients

        # Reference: Eichler, Koelmel, Sauer, J. of Comput. Chem., 18(4). 1997, 463-477.

        for cell_L, iat_cl_sys, iat_sys, r, iat_link in embed.linkatoms:

            # calculate the bond distance (r_bond) at the border

            xyz = embed[iat_sys].get_position() - embed[iat_cl_sys].get_position() + cell_L

            # calculate the bond lenght and the factor f

            rbond = np.sqrt(np.dot(xyz, xyz))

            f = r / rbond

            # receive the gradients for the link atom

            fL = f_cl[iat_link]

            # dot product fL, xyz

            # apply corrections for each direction

            for idir in xrange(3):

                # correct the atom in the system

                # correct the atom in the cluster

        return f_sys_lo

"""This module defines an ASE interface to Turbomole

http://www.turbomole.com/

"""

import os, sys, string

import numpy as np

class Turbomole(Calculator):

    def __init__(self, label='turbomole', calculate_energy="dscf", calculate_forces="grad", post_HF = False):

        self.label = label

        self.converged = False

    def execute(self, command):

        from subprocess import Popen, PIPE

        try:

            proc = Popen([command], shell=True, stderr=PIPE)

            error = proc.communicate()[1]

                raise OSError(error)

            print "TM command: ", command, "successfully executed"

        except OSError, e:

            sys.exit(1)

    def get_potential_energy(self, atoms):

        return self.forces.copy()

    def get_stress(self, atoms):

    def set_atoms(self, atoms):

        if self.atoms == atoms:

                energy_tmp = float(line.split()[1])

                if self.post_HF:

                    energy_tmp += float(line.split()[4])

        self.e_total = energy_tmp * Hartree

            print 'Gradients not found'

            raise RuntimeError("Please check TURBOMOLE gradients")

        iline += len(self.atoms) + 1

        nline -= 1

        for i in xrange(iline, nline):

            line = string.replace(lines[i],'D','E')

            #tmp=np.append(forces_tmp,np.array\

        # go over all pairs of atoms

        for iat_sys in xrange(nat_sys):

            xyz = positions_new[iat_sys]

            iat_cl = self.atom_map_sys_cl[iat_sys]

            if iat_cl > -1:

                self.atoms_cluster[iat_cl].set_position(xyz)

        for cell_L, iat_cl_sys, iat_sys, r, iat in self.linkatoms:

            # determine position of the link atom

            xyz *= r / np.sqrt(np.dot(xyz, xyz))

            xyz += xyz_cl

            # update xyz coordinates of the cluster

        f.close()

    atoms = Atoms(positions = atoms_pos, symbols = atom_symbols, pbc = False)

    #print c

            f.write('%20.14f  %20.14f  %20.14f      %2s \n' 

                    % (x/Bohr, y/Bohr, z/Bohr, s.lower()))

    f.write("$end\n")