<commit_before>// Copyright (c) 2013-2014 <NAME>, <NAME>
// All rights reserved.
// 
// Redistribution and use in source and binary forms, with or without modification, are permitted provided that 
// the following conditions are met:
// 
// 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the 
//    following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions 
//    and the following disclaimer in the documentation and/or other materials provided with the distribution.
// 
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED 
// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 
// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR 
// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

/** \file dft_ground_state.cpp
 *
 *  \brief Contains remaining implementation of sirius::DFT_ground_state class.
 */

#include "dft_ground_state.h"

namespace sirius {

double DFT_ground_state::energy_enuc()
{
    double enuc = 0.0;
    if (parameters_.unit_cell()->full_potential())
    {
        for (size_t ialoc = 0; ialoc < parameters_.unit_cell()->spl_num_atoms().local_size(); ialoc++)
        {
            int ia = parameters_.unit_cell()->spl_num_atoms((int)ialoc);
            int zn = parameters_.unit_cell()->atom(ia)->type()->zn();
            enuc -= 0.5 * zn * potential_->vh_el(ia) * y00;
        }
        parameters_.comm().allreduce(&enuc, 1);
    }
    
    return enuc;
}

double DFT_ground_state::core_eval_sum()
{
    double sum = 0.0;
    for (int ic = 0; ic < parameters_.unit_cell()->num_atom_symmetry_classes(); ic++)
    {
        sum += parameters_.unit_cell()->atom_symmetry_class(ic)->core_eval_sum() * 
               parameters_.unit_cell()->atom_symmetry_class(ic)->num_atoms();
    }
    return sum;
}

void DFT_ground_state::move_atoms(int istep)
{
    mdarray<double, 2> atom_force(3, parameters_.unit_cell()->num_atoms());
    forces(atom_force);
    if (verbosity_level >= 6 && parameters_.comm().rank() == 0)
    {
        printf("\n");
        printf("Atomic forces\n");
        for (int ia = 0; ia < parameters_.unit_cell()->num_atoms(); ia++)
        {
            printf("ia : %i, force : %12.6f %12.6f %12.6f\n", ia, atom_force(0, ia), atom_force(1, ia), atom_force(2, ia));
        }
    }

    for (int ia = 0; ia < parameters_.unit_cell()->num_atoms(); ia++)
    {
        vector3d<double> pos = parameters_.unit_cell()->atom(ia)->position();

        vector3d<double> forcef = parameters_.unit_cell()->get_fractional_coordinates(vector3d<double>(&atom_force(0, ia)));

        for (int x = 0; x < 3; x++) pos[x] += forcef[x];
        
        parameters_.unit_cell()->atom(ia)->set_position(pos);
    }
}

void DFT_ground_state::update()
{
    parameters_.update();
    potential_->update();
    kset_->update();
}

void DFT_ground_state::forces(mdarray<double, 2>& forces__)
{
    Force::total_force(parameters_, potential_, density_, kset_, forces__);
}

void DFT_ground_state::scf_loop(double potential_tol, double energy_tol, int num_dft_iter)
{
    Timer t("sirius::DFT_ground_state::scf_loop");
    
    density_->mixer_init();

    double eold = 0.0;
    double rms = 1.0;

    for (int iter = 0; iter < num_dft_iter; iter++)
    {
        Timer t1("sirius::DFT_ground_state::scf_loop|iteration");

        /* compute new potential */
        generate_effective_potential();

        //if (iter == 0) 
        //{
        //    potential_->mixer_init();
        //}
        //else
        //{
        //    potential_->mix();
        //}
        
        /* find new wave-functions */
        kset_->find_eigen_states(potential_, true);
        kset_->find_band_occupancies();

        /* generate new density from the occupied wave-functions */
        density_->generate(*kset_);

        if (use_symmetry_)
        {
            if (parameters_.esm_type() == full_potential_lapwlo || parameters_.esm_type() == full_potential_pwlo)
            {
                for (int j = 0; j < parameters_.num_mag_dims(); j++)
                    density_->magnetization(j)->fft_transform(-1);
            }

            symmetrize_density();

            if (parameters_.esm_type() == full_potential_lapwlo || parameters_.esm_type() == full_potential_pwlo)
            {
                density_->rho()->fft_transform(1);
                for (int j = 0; j < parameters_.num_mag_dims(); j++)
                    density_->magnetization(j)->fft_transform(1);
            }
        }
        
        rms = density_->mix();

        parameters_.comm().bcast(&rms, 1, 0);

        /* compute new total energy for a new density */
        double etot = total_energy();
        
        /* write some information */
        print_info();

        if (parameters_.comm().rank() == 0)
        {
            printf("iteration : %3i, density RMS %12.6f, energy difference : %12.6f\n", iter, rms, etot - eold);
        }
        
        if (std::abs(eold - etot) < energy_tol && rms < potential_tol) break;

        eold = etot;
    }
    
    parameters_.create_storage_file();
    potential_->save();
    density_->save();
}

void DFT_ground_state::relax_atom_positions()
{
    for (int i = 0; i < 5; i++)
    {
        scf_loop(1e-4, 1e-4, 100);
        move_atoms(i);
        update();<commit_msg>small fix<commit_after>// Copyright (c) 2013-2014 <NAME>, <NAME>
// All rights reserved.
// 
// Redistribution and use in source and binary forms, with or without modification, are permitted provided that 
// the following conditions are met:
// 
// 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the 
//    following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions 
//    and the following disclaimer in the documentation and/or other materials provided with the distribution.
// 
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED 
// WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 
// PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR 
// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

/** \file dft_ground_state.cpp
 *
 *  \brief Contains remaining implementation of sirius::DFT_ground_state class.
 */

#include "dft_ground_state.h"

namespace sirius {

double DFT_ground_state::energy_enuc()
{
    double enuc = 0.0;
    if (parameters_.unit_cell()->full_potential())
    {
        for (size_t ialoc = 0; ialoc < parameters_.unit_cell()->spl_num_atoms().local_size(); ialoc++)
        {
            int ia = parameters_.unit_cell()->spl_num_atoms((int)ialoc);
            int zn = parameters_.unit_cell()->atom(ia)->type()->zn();
            enuc -= 0.5 * zn * potential_->vh_el(ia) * y00;
        }
        parameters_.comm().allreduce(&enuc, 1);
    }
    
    return enuc;
}

double DFT_ground_state::core_eval_sum()
{
    double sum = 0.0;
    for (int ic = 0; ic < parameters_.unit_cell()->num_atom_symmetry_classes(); ic++)
    {
        sum += parameters_.unit_cell()->atom_symmetry_class(ic)->core_eval_sum() * 
               parameters_.unit_cell()->atom_symmetry_class(ic)->num_atoms();
    }
    return sum;
}

void DFT_ground_state::move_atoms(int istep)
{
    mdarray<double, 2> atom_force(3, parameters_.unit_cell()->num_atoms());
    forces(atom_force);
    if (verbosity_level >= 6 && parameters_.comm().rank() == 0)
    {
        printf("\n");
        printf("Atomic forces\n");
        for (int ia = 0; ia < parameters_.unit_cell()->num_atoms(); ia++)
        {
            printf("ia : %i, force : %12.6f %12.6f %12.6f\n", ia, atom_force(0, ia), atom_force(1, ia), atom_force(2, ia));
        }
    }

    for (int ia = 0; ia < parameters_.unit_cell()->num_atoms(); ia++)
    {
        vector3d<double> pos = parameters_.unit_cell()->atom(ia)->position();

        vector3d<double> forcef = parameters_.unit_cell()->get_fractional_coordinates(vector3d<double>(&atom_force(0, ia)));

        for (int x = 0; x < 3; x++) pos[x] += forcef[x];
        
        parameters_.unit_cell()->atom(ia)->set_position(pos);
    }
}

void DFT_ground_state::update()
{
    parameters_.update();
    potential_->update();
    kset_->update();
}

void DFT_ground_state::forces(mdarray<double, 2>& forces__)
{
    Force::total_force(parameters_, potential_, density_, kset_, forces__);
}

void DFT_ground_state::scf_loop(double potential_tol, double energy_tol, int num_dft_iter)
{
    Timer t("sirius::DFT_ground_state::scf_loop");
    
    density_->mixer_init();

    double eold = 0.0;

    for (int iter = 0; iter < num_dft_iter; iter++)
    {
        Timer t1("sirius::DFT_ground_state::scf_loop|iteration");

        /* compute new potential */
        generate_effective_potential();

        /* find new wave-functions */
        kset_->find_eigen_states(potential_, true);
        kset_->find_band_occupancies();

        /* generate new density from the occupied wave-functions */
        density_->generate(*kset_);

        if (use_symmetry_)
        {
            if (parameters_.esm_type() == full_potential_lapwlo || parameters_.esm_type() == full_potential_pwlo)
            {
                for (int j = 0; j < parameters_.num_mag_dims(); j++)
                    density_->magnetization(j)->fft_transform(-1);
            }

            symmetrize_density();

            if (parameters_.esm_type() == full_potential_lapwlo || parameters_.esm_type() == full_potential_pwlo)
            {
                density_->rho()->fft_transform(1);
                for (int j = 0; j < parameters_.num_mag_dims(); j++)
                    density_->magnetization(j)->fft_transform(1);
            }
        }
        
        double rms = density_->mix();

        parameters_.comm().bcast(&rms, 1, 0);

        /* compute new total energy for a new density */
        double etot = total_energy();
        
        /* write some information */
        print_info();

        if (parameters_.comm().rank() == 0)
        {
            printf("iteration : %3i, density RMS %12.6f, energy difference : %12.6f\n", iter, rms, etot - eold);
        }
        
        if (std::abs(eold - etot) < energy_tol && rms < potential_tol) break;

        eold = etot;
    }
    
    parameters_.create_storage_file();
    potential_->save();
    density_->save();
}

void DFT_ground_state::relax_atom_positions()
{
    for (int i = 0; i < 5; i++)
    {
        scf_loop(1e-4, 1e-4, 100);
        move_atoms(i);
        update();
{'fraction_non_alphanumeric': 0.13070717473395124, 'fraction_numerical': 0.012959148644009611, 'mean_word_length': 2.691162496040545, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 91, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '482488', 'n_tokens_mistral': 3775, 'n_tokens_neox': 3520, 'n_words': 926}