11#include " NCLibxc.h"
22#include < iostream>
33#include < iomanip>
4-
4+ namespace NCXC {
55// /////////////////////////////////////////////////////////////////////////////////
66// collinear limit test for GGA
77void NCLibxc::gga_collinear_test ()
@@ -640,4 +640,178 @@ void NCLibxc::gga_local_torque_test()
640640 << torque[3 *i+1 ] << " , "
641641 << torque[3 *i+2 ] << " )" << std::endl;
642642 }
643+ }
644+
645+ // xc local torque and xc potential from gga, prove well-defined torque and potential
646+ void NCLibxc::gga_deri_limit ()
647+ {
648+ // We now vary lambda from 1.0 down to 1e-20
649+ int xc_id = 106 ;
650+ std::complex <double > twoi (0.0 , 2.0 );
651+ std::complex <double > two (2.0 , 0.0 );
652+ for (int i = 0 ; i <= 20 ; ++i)
653+ {
654+ double lambda = std::pow (10.0 , -static_cast <double >(i));
655+
656+ // We only need one point
657+ std::vector<double > n (1 , 2.0 );
658+ std::vector<double > mx (1 , lambda);
659+ std::vector<double > my (1 , 0.1 *lambda);
660+ std::vector<double > mz (1 , 0.01 * lambda);
661+
662+ // All gradients/derivatives set to 0.01
663+ std::vector<double > gradx_n (1 , 0.01 ), grady_n (1 , 0.01 ), gradz_n (1 , 0.01 );
664+ std::vector<double > gradx_mx (1 , 0.01 ), grady_mx (1 , 0.01 ), gradz_mx (1 , 0.01 );
665+ std::vector<double > gradx_my (1 , 0.01 ), grady_my (1 , 0.01 ), gradz_my (1 , 0.01 );
666+ std::vector<double > gradx_mz (1 , 0.01 ), grady_mz (1 , 0.01 ), gradz_mz (1 , 0.01 );
667+
668+ std::vector<double > grad2xx_n (1 , 0.01 ), grad2yy_n (1 , 0.01 ), grad2zz_n (1 , 0.01 );
669+ std::vector<double > grad2xy_n (1 , 0.01 ), grad2yz_n (1 , 0.01 ), grad2xz_n (1 , 0.01 );
670+ std::vector<double > grad2xx_mx (1 , 0.01 ), grad2yy_mx (1 , 0.01 ), grad2zz_mx (1 , 0.01 );
671+ std::vector<double > grad2xy_mx (1 , 0.01 ), grad2yz_mx (1 , 0.01 ), grad2xz_mx (1 , 0.01 );
672+ std::vector<double > grad2xx_my (1 , 0.01 ), grad2yy_my (1 , 0.01 ), grad2zz_my (1 , 0.01 );
673+ std::vector<double > grad2xy_my (1 , 0.01 ), grad2yz_my (1 , 0.01 ), grad2xz_my (1 , 0.01 );
674+ std::vector<double > grad2xx_mz (1 , 0.01 ), grad2yy_mz (1 , 0.01 ), grad2zz_mz (1 , 0.01 );
675+ std::vector<double > grad2xy_mz (1 , 0.01 ), grad2yz_mz (1 , 0.01 ), grad2xz_mz (1 , 0.01 );
676+
677+ std::vector<double > grad3xxx_n (1 , 0.01 ), grad3xxy_n (1 , 0.01 ), grad3xxz_n (1 , 0.01 );
678+ std::vector<double > grad3xyy_n (1 , 0.01 ), grad3xyz_n (1 , 0.01 ), grad3xzz_n (1 , 0.01 );
679+ std::vector<double > grad3yyy_n (1 , 0.01 ), grad3yyz_n (1 , 0.01 ), grad3yzz_n (1 , 0.01 );
680+ std::vector<double > grad3zzz_n (1 , 0.01 );
681+
682+ std::vector<double > grad3xxx_mx (1 , 0.01 ), grad3xxy_mx (1 , 0.01 ), grad3xxz_mx (1 , 0.01 );
683+ std::vector<double > grad3xyy_mx (1 , 0.01 ), grad3xyz_mx (1 , 0.01 ), grad3xzz_mx (1 , 0.01 );
684+ std::vector<double > grad3yyy_mx (1 , 0.01 ), grad3yyz_mx (1 , 0.01 ), grad3yzz_mx (1 , 0.01 );
685+ std::vector<double > grad3zzz_mx (1 , 0.01 );
686+
687+ std::vector<double > grad3xxx_my (1 , 0.01 ), grad3xxy_my (1 , 0.01 ), grad3xxz_my (1 , 0.01 );
688+ std::vector<double > grad3xyy_my (1 , 0.01 ), grad3xyz_my (1 , 0.01 ), grad3xzz_my (1 , 0.01 );
689+ std::vector<double > grad3yyy_my (1 , 0.01 ), grad3yyz_my (1 , 0.01 ), grad3yzz_my (1 , 0.01 );
690+ std::vector<double > grad3zzz_my (1 , 0.01 );
691+
692+ std::vector<double > grad3xxx_mz (1 , 0.01 ), grad3xxy_mz (1 , 0.01 ), grad3xxz_mz (1 , 0.01 );
693+ std::vector<double > grad3xyy_mz (1 , 0.01 ), grad3xyz_mz (1 , 0.01 ), grad3xzz_mz (1 , 0.01 );
694+ std::vector<double > grad3yyy_mz (1 , 0.01 ), grad3yyz_mz (1 , 0.01 ), grad3yzz_mz (1 , 0.01 );
695+ std::vector<double > grad3zzz_mz (1 , 0.01 );
696+
697+ std::vector<double > grad4xxxx_n (1 , 0.01 ), grad4xxxy_n (1 , 0.01 ), grad4xxxz_n (1 , 0.01 );
698+ std::vector<double > grad4xxyy_n (1 , 0.01 ), grad4xxyz_n (1 , 0.01 ), grad4xxzz_n (1 , 0.01 );
699+ std::vector<double > grad4xyyy_n (1 , 0.01 ), grad4xyyz_n (1 , 0.01 ), grad4xyzz_n (1 , 0.01 );
700+ std::vector<double > grad4xzzz_n (1 , 0.01 ), grad4yyyy_n (1 , 0.01 ), grad4yyyz_n (1 , 0.01 );
701+ std::vector<double > grad4yyzz_n (1 , 0.01 ), grad4yzzz_n (1 , 0.01 ), grad4zzzz_n (1 , 0.01 );
702+
703+ std::vector<double > grad4xxxx_mx (1 , 0.01 ), grad4xxxy_mx (1 , 0.01 ), grad4xxxz_mx (1 , 0.01 );
704+ std::vector<double > grad4xxyy_mx (1 , 0.01 ), grad4xxyz_mx (1 , 0.01 ), grad4xxzz_mx (1 , 0.01 );
705+ std::vector<double > grad4xyyy_mx (1 , 0.01 ), grad4xyyz_mx (1 , 0.01 ), grad4xyzz_mx (1 , 0.01 );
706+ std::vector<double > grad4xzzz_mx (1 , 0.01 ), grad4yyyy_mx (1 , 0.01 ), grad4yyyz_mx (1 , 0.01 );
707+ std::vector<double > grad4yyzz_mx (1 , 0.01 ), grad4yzzz_mx (1 , 0.01 ), grad4zzzz_mx (1 , 0.01 );
708+
709+ std::vector<double > grad4xxxx_my (1 , 0.01 ), grad4xxxy_my (1 , 0.01 ), grad4xxxz_my (1 , 0.01 );
710+ std::vector<double > grad4xxyy_my (1 , 0.01 ), grad4xxyz_my (1 , 0.01 ), grad4xxzz_my (1 , 0.01 );
711+ std::vector<double > grad4xyyy_my (1 , 0.01 ), grad4xyyz_my (1 , 0.01 ), grad4xyzz_my (1 , 0.01 );
712+ std::vector<double > grad4xzzz_my (1 , 0.01 ), grad4yyyy_my (1 , 0.01 ), grad4yyyz_my (1 , 0.01 );
713+ std::vector<double > grad4yyzz_my (1 , 0.01 ), grad4yzzz_my (1 , 0.01 ), grad4zzzz_my (1 , 0.01 );
714+
715+ std::vector<double > grad4xxxx_mz (1 , 0.01 ), grad4xxxy_mz (1 , 0.01 ), grad4xxxz_mz (1 , 0.01 );
716+ std::vector<double > grad4xxyy_mz (1 , 0.01 ), grad4xxyz_mz (1 , 0.01 ), grad4xxzz_mz (1 , 0.01 );
717+ std::vector<double > grad4xyyy_mz (1 , 0.01 ), grad4xyyz_mz (1 , 0.01 ), grad4xyzz_mz (1 , 0.01 );
718+ std::vector<double > grad4xzzz_mz (1 , 0.01 ), grad4yyyy_mz (1 , 0.01 ), grad4yyyz_mz (1 , 0.01 );
719+ std::vector<double > grad4yyzz_mz (1 , 0.01 ), grad4yzzz_mz (1 , 0.01 ), grad4zzzz_mz (1 , 0.01 );
720+
721+ // Compute local torque
722+ auto torque = gga_torque (
723+ xc_id, n, mx, my, mz,
724+ gradx_n, grady_n, gradz_n,
725+ gradx_mx, grady_mx, gradz_mx,
726+ gradx_my, grady_my, gradz_my,
727+ gradx_mz, grady_mz, gradz_mz,
728+ grad2xx_n, grad2yy_n, grad2zz_n, grad2xy_n, grad2yz_n, grad2xz_n,
729+ grad2xx_mx, grad2yy_mx, grad2zz_mx, grad2xy_mx, grad2yz_mx, grad2xz_mx,
730+ grad2xx_my, grad2yy_my, grad2zz_my, grad2xy_my, grad2yz_my, grad2xz_my,
731+ grad2xx_mz, grad2yy_mz, grad2zz_mz, grad2xy_mz, grad2yz_mz, grad2xz_mz,
732+ grad3xxx_n, grad3xxy_n, grad3xxz_n, grad3xyy_n, grad3xyz_n, grad3xzz_n,
733+ grad3yyy_n, grad3yyz_n, grad3yzz_n, grad3zzz_n,
734+ grad3xxx_mx, grad3xxy_mx, grad3xxz_mx, grad3xyy_mx, grad3xyz_mx, grad3xzz_mx,
735+ grad3yyy_mx, grad3yyz_mx, grad3yzz_mx, grad3zzz_mx,
736+ grad3xxx_my, grad3xxy_my, grad3xxz_my, grad3xyy_my, grad3xyz_my, grad3xzz_my,
737+ grad3yyy_my, grad3yyz_my, grad3yzz_my, grad3zzz_my,
738+ grad3xxx_mz, grad3xxy_mz, grad3xxz_mz, grad3xyy_mz, grad3xyz_mz, grad3xzz_mz,
739+ grad3yyy_mz, grad3yyz_mz, grad3yzz_mz, grad3zzz_mz,
740+ grad4xxxx_n, grad4xxxy_n, grad4xxxz_n, grad4xxyy_n, grad4xxyz_n, grad4xxzz_n,
741+ grad4xyyy_n, grad4xyyz_n, grad4xyzz_n, grad4xzzz_n, grad4yyyy_n, grad4yyyz_n,
742+ grad4yyzz_n, grad4yzzz_n, grad4zzzz_n,
743+ grad4xxxx_mx, grad4xxxy_mx, grad4xxxz_mx, grad4xxyy_mx, grad4xxyz_mx, grad4xxzz_mx,
744+ grad4xyyy_mx, grad4xyyz_mx, grad4xyzz_mx, grad4xzzz_mx, grad4yyyy_mx, grad4yyyz_mx,
745+ grad4yyzz_mx, grad4yzzz_mx, grad4zzzz_mx,
746+ grad4xxxx_my, grad4xxxy_my, grad4xxxz_my, grad4xxyy_my, grad4xxyz_my, grad4xxzz_my,
747+ grad4xyyy_my, grad4xyyz_my, grad4xyzz_my, grad4xzzz_my, grad4yyyy_my, grad4yyyz_my,
748+ grad4yyzz_my, grad4yzzz_my, grad4zzzz_my,
749+ grad4xxxx_mz, grad4xxxy_mz, grad4xxxz_mz, grad4xxyy_mz, grad4xxyz_mz, grad4xxzz_mz,
750+ grad4xyyy_mz, grad4xyyz_mz, grad4xyzz_mz, grad4xzzz_mz, grad4yyyy_mz, grad4yyyz_mz,
751+ grad4yyzz_mz, grad4yzzz_mz, grad4zzzz_mz
752+ );
753+
754+ // Compute E and V using gga_mc
755+ auto [energy, potential] = gga_mc (
756+ xc_id, n, mx, my, mz,
757+ gradx_n, grady_n, gradz_n,
758+ gradx_mx, grady_mx, gradz_mx,
759+ gradx_my, grady_my, gradz_my,
760+ gradx_mz, grady_mz, gradz_mz,
761+ grad2xx_n, grad2yy_n, grad2zz_n, grad2xy_n, grad2yz_n, grad2xz_n,
762+ grad2xx_mx, grad2yy_mx, grad2zz_mx, grad2xy_mx, grad2yz_mx, grad2xz_mx,
763+ grad2xx_my, grad2yy_my, grad2zz_my, grad2xy_my, grad2yz_my, grad2xz_my,
764+ grad2xx_mz, grad2yy_mz, grad2zz_mz, grad2xy_mz, grad2yz_mz, grad2xz_mz,
765+ grad3xxx_n, grad3xxy_n, grad3xxz_n, grad3xyy_n, grad3xyz_n, grad3xzz_n,
766+ grad3yyy_n, grad3yyz_n, grad3yzz_n, grad3zzz_n,
767+ grad3xxx_mx, grad3xxy_mx, grad3xxz_mx, grad3xyy_mx, grad3xyz_mx, grad3xzz_mx,
768+ grad3yyy_mx, grad3yyz_mx, grad3yzz_mx, grad3zzz_mx,
769+ grad3xxx_my, grad3xxy_my, grad3xxz_my, grad3xyy_my, grad3xyz_my, grad3xzz_my,
770+ grad3yyy_my, grad3yyz_my, grad3yzz_my, grad3zzz_my,
771+ grad3xxx_mz, grad3xxy_mz, grad3xxz_mz, grad3xyy_mz, grad3xyz_mz, grad3xzz_mz,
772+ grad3yyy_mz, grad3yyz_mz, grad3yzz_mz, grad3zzz_mz,
773+ grad4xxxx_n, grad4xxxy_n, grad4xxxz_n, grad4xxyy_n, grad4xxyz_n, grad4xxzz_n,
774+ grad4xyyy_n, grad4xyyz_n, grad4xyzz_n, grad4xzzz_n, grad4yyyy_n, grad4yyyz_n,
775+ grad4yyzz_n, grad4yzzz_n, grad4zzzz_n,
776+ grad4xxxx_mx, grad4xxxy_mx, grad4xxxz_mx, grad4xxyy_mx, grad4xxyz_mx, grad4xxzz_mx,
777+ grad4xyyy_mx, grad4xyyz_mx, grad4xyzz_mx, grad4xzzz_mx, grad4yyyy_mx, grad4yyyz_mx,
778+ grad4yyzz_mx, grad4yzzz_mx, grad4zzzz_mx,
779+ grad4xxxx_my, grad4xxxy_my, grad4xxxz_my, grad4xxyy_my, grad4xxyz_my, grad4xxzz_my,
780+ grad4xyyy_my, grad4xyyz_my, grad4xyzz_my, grad4xzzz_my, grad4yyyy_my, grad4yyyz_my,
781+ grad4yyzz_my, grad4yzzz_my, grad4zzzz_my,
782+ grad4xxxx_mz, grad4xxxy_mz, grad4xxxz_mz, grad4xxyy_mz, grad4xxyz_mz, grad4xxzz_mz,
783+ grad4xyyy_mz, grad4xyyz_mz, grad4xyzz_mz, grad4xzzz_mz, grad4yyyy_mz, grad4yyyz_mz,
784+ grad4yyzz_mz, grad4yzzz_mz, grad4zzzz_mz
785+ );
786+
787+ // Print results
788+ std::cout << std::fixed << std::setprecision (30 );
789+ std::cout << " For lambda=" << lambda << " :\n " ;
790+
791+ // Torque
792+ std::cout << " Torque: ("
793+ << torque[0 ] << " , "
794+ << torque[1 ] << " , "
795+ << torque[2 ] << " )" << std::endl;
796+
797+ // Energy
798+ std::cout << " E = " << n[0 ]*energy[0 ] << std::endl;
799+
800+
801+ // Potential V: 2x2 matrix
802+ const auto & mat = potential[0 ];
803+ double V0 = std::real ((mat[0 ][0 ]+mat[1 ][1 ])/two);
804+ double V1 = std::real ((mat[0 ][1 ]+mat[1 ][0 ])/two);
805+ double V2 = std::real ((mat[1 ][0 ]-mat[0 ][1 ])/twoi);
806+ double V3 = std::real ((mat[0 ][0 ] - mat[1 ][1 ]) / two);
807+ std::cout << " DE/Dn =" << V0
808+ << " , DE/Dmx = " << V1
809+ << " , DE/Dmy = " << V2
810+ << " , DE/Dmz = " << V3 << std::endl;
811+
812+
813+ std::cout << std::endl;
814+ }
815+ }
816+
643817}
0 commit comments