
% ------------------------------------------------------------------------
% Function: minfunc_svpn ( Input )
% ------------------------------------------------------------------------
% Aim:
% The total energy of a dislocation according the SVPN model
% ------------------------------------------------------------------------
% Ref:
% 1/ Bulatov, Vasily V., and Efthimios Kaxiras. Physical Review Letters 78.22 (1997): 4221. 
% 2/ Lu, Gang. Handbook of Materials Modeling. Springer Netherlands, 2005. 793-811.
% ------------------------------------------------------------------------

function MIN = minfunc_svpn ( Input )

    %   defining global variables

    global Dx Wx Dz Wz Ax Az;
    global BurVect;
    global shear_modulus Kx Kz;
    global dx dz;
    global Nmis;
    global dax;
    global Inpas;
    global StressEnergy GradientEnergy NonlocalEnergy;
    global str_coef grad_coef nonloc_coef;
    
    %   outer cut-off radius

    Rcut=50;
    
    %% >>>>>> enter into main loop <<<<<<

    nx=size(Input,1);
    
    %   To define some Matrix
    
    MIN=zeros(nx,1);
    
    for i=1:nx
    
        if Nmis>=1
            
            % % x direction
            
            for j=1:(Nmis-1)
                Ax(j)=Input(i,j);
            end  
            for j=1:(Nmis-1)
                Dx(j)=Input(i,Nmis-1+j);
            end         
            for j=1:Nmis
                Wx(j)=Input(i,2*Nmis-2+j);
            end
            
            % % z direction
    
            for j=1:(Nmis-1)
                Az(j)=Input(i,3*Nmis-2+j);
            end
            for j=1:Nmis-1
                Dz(j)=Input(i,4*Nmis-3+j);
            end           
            for j=1:Nmis
                Wz(j)=Input(i,5*Nmis-4+j);
            end            
            
        elseif Nmis==0

            %   some parameters

            Ax=1/2;
            Az=sqrt(3)/6;

            Dx=Input(i,1);
            Wx=Input(i,2);
            Dz=Input(i,3);
            Wz=Input(i,4);
            
        elseif Nmis==-1

            Ax=Input(i,1);
            Dx=Input(i,2);
            Wx=Input(i,3);
            dx=Input(i,4);
            
            Az=Input(i,5);
            Dz=Input(i,6);
            Wz=Input(i,7);
            dz=Input(i,8);
            
        end

        %% elastic energy (1e-11*J/m)

        [n,nn]=meshgrid(((-Rcut+dax):(Rcut+dax)),((-Rcut+dax):(Rcut+dax)));

        Chi=3/2*Inpas^2+1/2*((n-1)*Inpas-(nn-1)*Inpas).^2.*log(abs((n-1)*Inpas-(nn-1)*Inpas))...
                      +1/2*(n*Inpas-nn*Inpas).^2.*log(abs(n*Inpas-nn*Inpas))...
                      -1/2*(n*Inpas-(nn-1)*Inpas).^2.*log(abs(n*Inpas-(nn-1)*Inpas))...
                      -1/2*((n-1)*Inpas-nn*Inpas).^2.*log(abs((n-1)*Inpas-nn*Inpas));

        Chi(1,1)=3/2*Inpas^2-Inpas^2*log(Inpas);
        for m=2:(2*Rcut+1)
            Chi(m,m)=3/2*Inpas^2-Inpas^2*log(Inpas);
            Chi(m,m-1)=3/2*Inpas^2+Inpas^2*log(Inpas)-1/2*(2*Inpas)^2*log(2*Inpas);
            Chi(m-1,m)=Chi(m,m-1);
        end

        % Eel=1/(4*pi)*Kx*BurVect^2*log(Rcut*BurVect)+sum(sum(Chi.*(1/(4*pi)*Kx*(misfit_x(n*Inpas)-misfit_x((n-1)*Inpas))/Inpas.*(misfit_x(nn*Inpas)-misfit_x((nn-1)*Inpas))/Inpas...
        %                                                          +1/(4*pi)*Kz*(misfit_z(n*Inpas)-misfit_z((n-1)*Inpas))/Inpas.*(misfit_z(nn*Inpas)-misfit_z((nn-1)*Inpas))/Inpas)));

        Eel=sum(sum(Chi.*(1/(4*pi)*Kx*(misfit_x(n*Inpas)-misfit_x((n-1)*Inpas))/Inpas.*(misfit_x(nn*Inpas)-misfit_x((nn-1)*Inpas))/Inpas...
                         +1/(4*pi)*Kz*(misfit_z(n*Inpas)-misfit_z((n-1)*Inpas))/Inpas.*(misfit_z(nn*Inpas)-misfit_z((nn-1)*Inpas))/Inpas)));        
        
        %% misfit energy (1e-10*J/m)

        mm=-(Rcut*10+dax):(Rcut*10+dax);
        EA=sum((Gamma(mm*Inpas)*Inpas));

        %% stress interaction term (1e-10*J/m and unit of str_coef: MPa)
    
        if StressEnergy==0
    
            Estr=0.0000;
            
        elseif StressEnergy==1
            
            % Estr=-1E-4*str_coef*sum(((mm*Inpas).^2-((mm-1)*Inpas).^2)/2.*(misfit_x(mm*Inpas)-misfit_x((mm-1)*Inpas))/Inpas);
            Estr=-1E-4*str_coef*sum((misfit_x(mm*Inpas)+misfit_x((mm+1)*Inpas))*Inpas)/2;
        end
        
        %% gradient energy (1e-10*J/m and unit of grad_coef: shear_modulus*Burvect)
    
        if GradientEnergy==0
            
            Egradient=0.0000;
            
        elseif GradientEnergy==1
            
            Egradient=1E-1*grad_coef*shear_modulus*BurVect*Inpas*(sum(((misfit_x((mm+1)*Inpas)-misfit_x((mm-1)*Inpas))/(2*Inpas)).^2)...
                                                                 +sum(((misfit_z((mm+1)*Inpas)-misfit_z((mm-1)*Inpas))/(2*Inpas)).^2));
                                  
        end
                              
        %% simplified nonlocal energy (1e-10*J/m and unit of nonlocal_coef: shear_modulus/Burvect)
        
        if NonlocalEnergy==0
            
            Enonlocal=0.0000;
        
        elseif NonlocalEnergy==1
            
            Enonlocal=1E-1*nonloc_coef*shear_modulus/BurVect*(sum(misfit_x(mm*Inpas).*(misfit_x(mm*Inpas)-(misfit_x((mm+1)*Inpas)+misfit_x((mm-1)*Inpas))/2)*Inpas)...
                                                             +sum(misfit_z(mm*Inpas).*(misfit_z(mm*Inpas)-(misfit_z((mm+1)*Inpas)+misfit_z((mm-1)*Inpas))/2)*Inpas));
        
        end
                          
        %% total energy(unit: 1e-10*J/m) of a dislocation

        MIN(i,1)=1E-1*Eel + EA + Estr + Egradient + Enonlocal;

    end
    
end
