<commit_before>/*
 *            Copyright 2009-2012 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

// Overload of uBLAS prod function with MKL/GSL implementations
#include <votca/ctp/votca_ctp_config.h>

#include <votca/ctp/gwbse.h>

#include <boost/format.hpp>
#include <boost/filesystem.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <votca/ctp/aomatrix.h>
#include <votca/ctp/threecenters.h>
// #include <votca/ctp/logger.h>
#include <votca/ctp/qmpackagefactory.h>
#include <boost/math/constants/constants.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <votca/tools/linalg.h>

using boost::format;
using namespace boost::filesystem;

namespace votca {
    namespace ctp {
        namespace ub = boost::numeric::ublas;

        // +++++++++++++++++++++++++++++ //
        // MBPT MEMBER FUNCTIONS         //
        // +++++++++++++++++++++++++++++ //

        
        

        void GWBSE::FullQPHamiltonian(){
            
            // constructing full QP Hamiltonian, storage in vxc
            _vxc = -_vxc + _sigma_x + _sigma_c;
            // diagonal elements are given by _qp_energies
            for (int _m = 0; _m < _vxc.size1(); _m++ ){
              _vxc( _m,_m ) = _qp_energies( _m + _qpmin );
            }

            
            // sigma matrices can be freed
            _sigma_x.resize(0,0);
            _sigma_c.resize(0,0);
            
            
            if ( _do_qp_diag ){
            /* diagonalize, since _vxc will be needed in BSE, and GSL
             * destroys the input array, we need to make a local copy first
             */
            
            // get eigenvalues and eigenvectors of this matrix
            ub::matrix<double> _temp = _vxc;
            _qp_diag_energies.resize(_temp.size1());
            _qp_diag_coefficients.resize(_temp.size1(), _temp.size1());
            linalg_eigenvalues(_temp, _qp_diag_energies, _qp_diag_coefficients);

            // TODO storage -> orbitals

            }
            
        }
        
        
        void GWBSE::sigma_c_setup(const TCMatrix& _Mmn, const ub::vector<double>& _edft){
            
            // iterative refinement of qp energies
            int _max_iter = 5;
            int _levelsum = _Mmn[0].size2(); // total number of bands
            int _gwsize  = _Mmn[0].size1(); // size of the GW basis
            const double pi = boost::math::constants::pi<double>();
            

            // initial _qp_energies are dft energies
            _qp_energies = _edft; // RANGES!
            double _DFTgap = _qp_energies( _homo +1 ) - _qp_energies( _homo  );
     

	    // only diagonal elements except for in final iteration
            for ( int _i_iter = 0 ; _i_iter < _max_iter-1 ; _i_iter++ ){
                
	      // initialize sigma_c to zero at the beginning of each iteration
	      _sigma_c = ub::zero_matrix<double>(_qptotal,_qptotal);

	      // loop over all GW levels
              #pragma omp parallel for
	      for (int _gw_level = 0; _gw_level < _qptotal ; _gw_level++ ){
                  
                const ub::matrix<double>& Mmn = _Mmn[ _gw_level + _qpmin ];
              
		// loop over all functions in GW basis
		for ( int _i_gw = 0; _i_gw < _gwsize ; _i_gw++ ){
                    
		  // loop over all bands
		  for ( int _i = 0; _i < _levelsum ; _i++ ){
                    
		    double occ = 1.0;
		    if ( _i > _homo ) occ = -1.0; // sign for empty levels
                                                    
		    // energy denominator
		    double _denom = _qp_energies( _gw_level + _qpmin  ) - _qp_energies( _i ) + occ*_ppm_freq( _i_gw );
                    
		    double _stab = 1.0;
		    if ( std::abs(_denom) < 0.5 ) {
		      _stab = 0.5 * ( 1.0 - cos(2.0 * pi * std::abs(_denom) ) );
		    }
                            
		    double _factor = _ppm_weight( _i_gw ) * _ppm_freq( _i_gw) * _stab/_denom; // contains conversion factor 2!
		    
		    // sigma_c diagonal elements
		    _sigma_c( _gw_level , _gw_level ) += _factor * Mmn(_i_gw, _i) *  Mmn(_i_gw, _i);  
                            
		  }// bands
                        
		}// GW functions
		
		// update _qp_energies
		_qp_energies( _gw_level  + _qpmin) = _edft( _gw_level + _qpmin ) + _sigma_x(_gw_level, _gw_level) + _sigma_c(_gw_level,_gw_level) - _vxc(_gw_level  ,_gw_level );
                    
	      }// all bands
                
            } // iterations
            
            // check HOMO-LUMO gap shift
            // cout << " QP HOMO : " << _qp_energies( _homo  ) << endl;
            // cout << " QP LUMO : " << _qp_energies( _homo +1 ) << endl;
            double _QPgap = _qp_energies( _homo +1 ) - _qp_energies( _homo  );
            // cout << " QP Gap  : " << _QPgap << endl;
            double _shift_new = _QPgap - _DFTgap;
            // cout << " shift new " << _shift_new << endl;
            if ( std::abs( (_shift_new - _shift)*13.605698066 ) > 0.01 ) {
                _shift = _shift_new;
            } else {
                _shift_converged = true;
            }

	    if ( ! _iterate_shift  ) _shift_converged = true;
            
            // only if _shift is converged
            if ( _shift_converged ){
	    // in final step, also calc offdiagonal elements
	    // initialize sigma_c to zero at the beginning
	    _sigma_c = ub::zero_matrix<double>(_qptotal,_qptotal);

            // loop over col  GW levels
            #pragma omp parallel for
	    for (int _gw_level = 0; _gw_level < _qptotal ; _gw_level++ ){
              
                const ub::matrix<double>& Mmn =  _Mmn[ _gw_level + _qpmin ];
<commit_msg>QP shift output<commit_after>/*
 *            Copyright 2009-2012 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

// Overload of uBLAS prod function with MKL/GSL implementations
#include <votca/ctp/votca_ctp_config.h>

#include <votca/ctp/gwbse.h>

#include <boost/format.hpp>
#include <boost/filesystem.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <votca/ctp/aomatrix.h>
#include <votca/ctp/threecenters.h>
// #include <votca/ctp/logger.h>
#include <votca/ctp/qmpackagefactory.h>
#include <boost/math/constants/constants.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <votca/tools/linalg.h>

using boost::format;
using namespace boost::filesystem;

namespace votca {
    namespace ctp {
        namespace ub = boost::numeric::ublas;

        // +++++++++++++++++++++++++++++ //
        // MBPT MEMBER FUNCTIONS         //
        // +++++++++++++++++++++++++++++ //

        
        

        void GWBSE::FullQPHamiltonian(){
            
            // constructing full QP Hamiltonian, storage in vxc
            _vxc = -_vxc + _sigma_x + _sigma_c;
            // diagonal elements are given by _qp_energies
            for (int _m = 0; _m < _vxc.size1(); _m++ ){
              _vxc( _m,_m ) = _qp_energies( _m + _qpmin );
            }

            
            // sigma matrices can be freed
            _sigma_x.resize(0,0);
            _sigma_c.resize(0,0);
            
            
            if ( _do_qp_diag ){
            /* diagonalize, since _vxc will be needed in BSE, and GSL
             * destroys the input array, we need to make a local copy first
             */
            
            // get eigenvalues and eigenvectors of this matrix
            ub::matrix<double> _temp = _vxc;
            _qp_diag_energies.resize(_temp.size1());
            _qp_diag_coefficients.resize(_temp.size1(), _temp.size1());
            linalg_eigenvalues(_temp, _qp_diag_energies, _qp_diag_coefficients);

            // TODO storage -> orbitals

            }
            
        }
        
        
        void GWBSE::sigma_c_setup(const TCMatrix& _Mmn, const ub::vector<double>& _edft){
            
            // iterative refinement of qp energies
            int _max_iter = 5;
            int _levelsum = _Mmn[0].size2(); // total number of bands
            int _gwsize  = _Mmn[0].size1(); // size of the GW basis
            const double pi = boost::math::constants::pi<double>();
            

            // initial _qp_energies are dft energies
            _qp_energies = _edft; // RANGES!
            double _DFTgap = _qp_energies( _homo +1 ) - _qp_energies( _homo  );
     

	    // only diagonal elements except for in final iteration
            for ( int _i_iter = 0 ; _i_iter < _max_iter-1 ; _i_iter++ ){
                
	      // initialize sigma_c to zero at the beginning of each iteration
	      _sigma_c = ub::zero_matrix<double>(_qptotal,_qptotal);

	      // loop over all GW levels
              #pragma omp parallel for
	      for (int _gw_level = 0; _gw_level < _qptotal ; _gw_level++ ){
                  
                const ub::matrix<double>& Mmn = _Mmn[ _gw_level + _qpmin ];
              
		// loop over all functions in GW basis
		for ( int _i_gw = 0; _i_gw < _gwsize ; _i_gw++ ){
                    
		  // loop over all bands
		  for ( int _i = 0; _i < _levelsum ; _i++ ){
                    
		    double occ = 1.0;
		    if ( _i > _homo ) occ = -1.0; // sign for empty levels
                                                    
		    // energy denominator
		    double _denom = _qp_energies( _gw_level + _qpmin  ) - _qp_energies( _i ) + occ*_ppm_freq( _i_gw );
                    
		    double _stab = 1.0;
		    if ( std::abs(_denom) < 0.5 ) {
		      _stab = 0.5 * ( 1.0 - cos(2.0 * pi * std::abs(_denom) ) );
		    }
                            
		    double _factor = _ppm_weight( _i_gw ) * _ppm_freq( _i_gw) * _stab/_denom; // contains conversion factor 2!
		    
		    // sigma_c diagonal elements
		    _sigma_c( _gw_level , _gw_level ) += _factor * Mmn(_i_gw, _i) *  Mmn(_i_gw, _i);  
                            
		  }// bands
                        
		}// GW functions
		
		// update _qp_energies
		_qp_energies( _gw_level  + _qpmin) = _edft( _gw_level + _qpmin ) + _sigma_x(_gw_level, _gw_level) + _sigma_c(_gw_level,_gw_level) - _vxc(_gw_level  ,_gw_level );
                    
	      }// all bands
                
            } // iterations
            
            // check HOMO-LUMO gap shift
            // cout << " QP HOMO : " << _qp_energies( _homo  ) << endl;
            // cout << " QP LUMO : " << _qp_energies( _homo +1 ) << endl;
            double _QPgap = _qp_energies( _homo +1 ) - _qp_energies( _homo  );
            // cout << " QP Gap  : " << _QPgap << endl;
            double _shift_new = _QPgap - _DFTgap;
            cout << " shift new " << _shift_new << endl;
            if ( std::abs( (_shift_new - _shift)*13.605698066 ) > 0.01 ) {
                _shift = _shift_new;
            } else {
                _shift_converged = true;
            }

	    if ( ! _iterate_shift  ) _shift_converged = true;
            
            // only if _shift is converged
            if ( _shift_converged ){
	    // in final step, also calc offdiagonal elements
	    // initialize sigma_c to zero at the beginning
	    _sigma_c = ub::zero_matrix<double>(_qptotal,_qptotal);

            // loop over col  GW levels
            #pragma omp parallel for
	    for (int _gw_level = 0; _gw_level < _qptotal ; _gw_level++ ){
              
                const ub::matrix<double>& Mmn =  _Mmn[ _gw_level + _qpmin ];

{'fraction_non_alphanumeric': 0.11067447269227684, 'fraction_numerical': 0.010265458336674953, 'mean_word_length': 1.693304535637149, 'pattern_counts': {'":': 0, '<': 103, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '82971', 'n_tokens_mistral': 4227, 'n_tokens_neox': 3929, 'n_words': 1022}