Phonon.h

/*-------------------------------------------------------------------
Copyright 2014 Ravishankar Sundararaman

This file is part of JDFTx.

JDFTx is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

JDFTx is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with JDFTx.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------*/

#ifndef JDFTX_PHONON_PHONON_H
#define JDFTX_PHONON_PHONON_H

#include <electronic/Everything.h>
#include <electronic/ColumnBundle.h>
#include <core/LatticeUtils.h>

class BlockRotationMatrix
{       int nBlocks, blockSize;
        std::vector<int> colOffset; 
        std::vector<matrix> rots; 
public:
        void init(int nBlocks, int blockSize); 
        void set(int rowBlock, int colBlock, const matrix& rot); 
        void allReduce(); 
        
        matrix transform(const matrix& in) const; 
};

struct PhononEverything: public Everything
{       class Phonon& phonon;
        PhononEverything(class Phonon& phonon);
};

class Phonon
{
public:
        std::vector<std::pair<string,string> > input; //input file contents
        
        vector3<int> sup; //phonon supercell 
        double dr; //perturbation amplitude in Cartesian coordinates
        double T; //temperature for free energy estimation
        double Fcut; //fillings cutoff for optimizing number of bands
        
        Phonon();
        void setup(bool printDefaults); //setup unit cell and basis modes for perturbations
        void dump(); //main calculations (sequence of supercell calculations) as well as output
        
private:
        PhononEverything e; //data for original unit cell
        PhononEverything eSupTemplate; //uninitialized version of eSup, with various flags later used to create eSup for each mode
        std::shared_ptr<PhononEverything> eSup; //supercell data for current perturbation

        int nSpins, nSpinor; //number of explicit spins and spinor length
        int nBandsOpt; //optimized number of bands, accounting for Fcut
        int prodSup; //number of unit cells in supercell
        std::vector< matrix3<> > symSupCart; //Cartesian symmetry rotation matrices for unperturbed supercell
        std::vector< std::vector<BlockRotationMatrix> > stateRot; //Unitary rotation of states involved in gamma-point Hsub for each supercell symmetry operation
        
        //Basis for phonon modes (not reduced by symmetries):
        struct Mode
        {       int sp, at; 
                vector3<> dir; 
        };
        std::vector<Mode> modes;
        std::vector<IonicGradient> dgrad; 
        std::vector< std::vector<matrix> > dHsub; 
        
        //Minimal basis of perturbations for supercell calculations:
        struct Perturbation : public Mode
        {       double weight; 
        };
        std::vector<Perturbation> perturbations;
        
        //Run supercell calculation for specified perturbation
        void processPerturbation(const Perturbation& pert);
        
        //set unperturbed state of supercell from unit cell
        //optionally get the subspace Hamiltonian at supercell Gamma point (for all bands, not just those in nBandsOpt)
        void setSupState(std::vector<matrix>* Hsub=0); 
        
        //Mapping between unit cell and current supercell k-points (generated by processPerturbation and used by setSupState)
        struct StateMapEntry : public Supercell::KmeshTransform //contain source k-point rotation here
        {       int qSup; //state index for supercell
                int nqPrev; //number of previous unit cell k-points that point to this supercell
                std::shared_ptr<class ColumnBundleTransform> transform; //wavefunction map
        };
        std::vector<StateMapEntry> stateMap; //map from unit cell k-points to supercell k-points
        
        //Utilities for mapping 3-dimensional and flattened cell indices:
        vector3<int> getCell(int unit) const;
        int getUnit(const vector3<int>& cell) const;
};

#endif //JDFTX_PHONON_PHONON_H