@@ -1020,282 +1020,6 @@ void UnitCell::setup(const std::string& latname_in,
10201020 return ;
10211021}
10221022
1023- void UnitCell::remake_cell () {
1024- ModuleBase::TITLE (" UnitCell" , " rmake_cell" );
1025-
1026- // The idea is as follows: for each type of lattice, first calculate
1027- // from current latvec the lattice parameters, then use the parameters
1028- // to reconstruct latvec
1029-
1030- if (latName == " none" ) {
1031- ModuleBase::WARNING_QUIT (
1032- " UnitCell" ,
1033- " to use fixed_ibrav, latname must be provided" );
1034- } else if (latName == " sc" ) // ibrav = 1
1035- {
1036- double celldm = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
1037- + pow (latvec.e13 , 2 ));
1038-
1039- latvec.Zero ();
1040- latvec.e11 = latvec.e22 = latvec.e33 = celldm;
1041- } else if (latName == " fcc" ) // ibrav = 2
1042- {
1043- double celldm = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
1044- + pow (latvec.e13 , 2 ))
1045- / std::sqrt (2.0 );
1046-
1047- latvec.e11 = -celldm;
1048- latvec.e12 = 0.0 ;
1049- latvec.e13 = celldm;
1050- latvec.e21 = 0.0 ;
1051- latvec.e22 = celldm;
1052- latvec.e23 = celldm;
1053- latvec.e31 = -celldm;
1054- latvec.e32 = celldm;
1055- latvec.e33 = 0.0 ;
1056- } else if (latName == " bcc" ) // ibrav = 3
1057- {
1058- double celldm = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
1059- + pow (latvec.e13 , 2 ))
1060- / std::sqrt (3.0 );
1061-
1062- latvec.e11 = celldm;
1063- latvec.e12 = celldm;
1064- latvec.e13 = celldm;
1065- latvec.e21 = -celldm;
1066- latvec.e22 = celldm;
1067- latvec.e23 = celldm;
1068- latvec.e31 = -celldm;
1069- latvec.e32 = -celldm;
1070- latvec.e33 = celldm;
1071- } else if (latName == " hexagonal" ) // ibrav = 4
1072- {
1073- double celldm1 = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
1074- + pow (latvec.e13 , 2 ));
1075- double celldm3 = std::sqrt (pow (latvec.e31 , 2 ) + pow (latvec.e32 , 2 )
1076- + pow (latvec.e33 , 2 ));
1077- double e22 = sqrt (3.0 ) / 2.0 ;
1078-
1079- latvec.e11 = celldm1;
1080- latvec.e12 = 0.0 ;
1081- latvec.e13 = 0.0 ;
1082- latvec.e21 = -0.5 * celldm1;
1083- latvec.e22 = celldm1 * e22 ;
1084- latvec.e23 = 0.0 ;
1085- latvec.e31 = 0.0 ;
1086- latvec.e32 = 0.0 ;
1087- latvec.e33 = celldm3;
1088- } else if (latName == " trigonal" ) // ibrav = 5
1089- {
1090- double celldm1 = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
1091- + pow (latvec.e13 , 2 ));
1092- double celldm2 = std::sqrt (pow (latvec.e21 , 2 ) + pow (latvec.e22 , 2 )
1093- + pow (latvec.e23 , 2 ));
1094- double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22
1095- + latvec.e13 * latvec.e23 );
1096- double cos12 = celldm12 / celldm1 / celldm2;
1097-
1098- if (cos12 <= -0.5 || cos12 >= 1.0 ) {
1099- ModuleBase::WARNING_QUIT (" unitcell" , " wrong cos12!" );
1100- }
1101- double t1 = sqrt (1.0 + 2.0 * cos12);
1102- double t2 = sqrt (1.0 - cos12);
1103-
1104- double e11 = celldm1 * t2 / sqrt (2.0 );
1105- double e12 = -celldm1 * t2 / sqrt (6.0 );
1106- double e13 = celldm1 * t1 / sqrt (3.0 );
1107- double e22 = celldm1 * sqrt (2.0 ) * t2 / sqrt (3.0 );
1108-
1109- latvec.e11 = e11 ;
1110- latvec.e12 = e12 ;
1111- latvec.e13 = e13 ;
1112- latvec.e21 = 0.0 ;
1113- latvec.e22 = e22 ;
1114- latvec.e23 = e13 ;
1115- latvec.e31 = -e11 ;
1116- latvec.e32 = e12 ;
1117- latvec.e33 = e13 ;
1118- } else if (latName == " st" ) // ibrav = 6
1119- {
1120- double celldm1 = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
1121- + pow (latvec.e13 , 2 ));
1122- double celldm3 = std::sqrt (pow (latvec.e31 , 2 ) + pow (latvec.e32 , 2 )
1123- + pow (latvec.e33 , 2 ));
1124- latvec.e11 = celldm1;
1125- latvec.e12 = 0.0 ;
1126- latvec.e13 = 0.0 ;
1127- latvec.e21 = 0.0 ;
1128- latvec.e22 = celldm1;
1129- latvec.e23 = 0.0 ;
1130- latvec.e31 = 0.0 ;
1131- latvec.e32 = 0.0 ;
1132- latvec.e33 = celldm3;
1133- } else if (latName == " bct" ) // ibrav = 7
1134- {
1135- double celldm1 = std::abs (latvec.e11 );
1136- double celldm2 = std::abs (latvec.e13 );
1137-
1138- latvec.e11 = celldm1;
1139- latvec.e12 = -celldm1;
1140- latvec.e13 = celldm2;
1141- latvec.e21 = celldm1;
1142- latvec.e22 = celldm1;
1143- latvec.e23 = celldm2;
1144- latvec.e31 = -celldm1;
1145- latvec.e32 = -celldm1;
1146- latvec.e33 = celldm2;
1147- } else if (latName == " so" ) // ibrav = 8
1148- {
1149- double celldm1 = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
1150- + pow (latvec.e13 , 2 ));
1151- double celldm2 = std::sqrt (pow (latvec.e21 , 2 ) + pow (latvec.e22 , 2 )
1152- + pow (latvec.e23 , 2 ));
1153- double celldm3 = std::sqrt (pow (latvec.e31 , 2 ) + pow (latvec.e32 , 2 )
1154- + pow (latvec.e33 , 2 ));
1155-
1156- latvec.e11 = celldm1;
1157- latvec.e12 = 0.0 ;
1158- latvec.e13 = 0.0 ;
1159- latvec.e21 = 0.0 ;
1160- latvec.e22 = celldm2;
1161- latvec.e23 = 0.0 ;
1162- latvec.e31 = 0.0 ;
1163- latvec.e32 = 0.0 ;
1164- latvec.e33 = celldm3;
1165- } else if (latName == " baco" ) // ibrav = 9
1166- {
1167- double celldm1 = std::abs (latvec.e11 );
1168- double celldm2 = std::abs (latvec.e22 );
1169- double celldm3 = std::abs (latvec.e33 );
1170-
1171- latvec.e11 = celldm1;
1172- latvec.e12 = celldm2;
1173- latvec.e13 = 0.0 ;
1174- latvec.e21 = -celldm1;
1175- latvec.e22 = celldm2;
1176- latvec.e23 = 0.0 ;
1177- latvec.e31 = 0.0 ;
1178- latvec.e32 = 0.0 ;
1179- latvec.e33 = celldm3;
1180- } else if (latName == " fco" ) // ibrav = 10
1181- {
1182- double celldm1 = std::abs (latvec.e11 );
1183- double celldm2 = std::abs (latvec.e22 );
1184- double celldm3 = std::abs (latvec.e33 );
1185-
1186- latvec.e11 = celldm1;
1187- latvec.e12 = 0.0 ;
1188- latvec.e13 = celldm3;
1189- latvec.e21 = celldm1;
1190- latvec.e22 = celldm2;
1191- latvec.e23 = 0.0 ;
1192- latvec.e31 = 0.0 ;
1193- latvec.e32 = celldm2;
1194- latvec.e33 = celldm3;
1195- } else if (latName == " bco" ) // ibrav = 11
1196- {
1197- double celldm1 = std::abs (latvec.e11 );
1198- double celldm2 = std::abs (latvec.e12 );
1199- double celldm3 = std::abs (latvec.e13 );
1200-
1201- latvec.e11 = celldm1;
1202- latvec.e12 = celldm2;
1203- latvec.e13 = celldm3;
1204- latvec.e21 = -celldm1;
1205- latvec.e22 = celldm2;
1206- latvec.e23 = celldm3;
1207- latvec.e31 = -celldm1;
1208- latvec.e32 = -celldm2;
1209- latvec.e33 = celldm3;
1210- } else if (latName == " sm" ) // ibrav = 12
1211- {
1212- double celldm1 = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
1213- + pow (latvec.e13 , 2 ));
1214- double celldm2 = std::sqrt (pow (latvec.e21 , 2 ) + pow (latvec.e22 , 2 )
1215- + pow (latvec.e23 , 2 ));
1216- double celldm3 = std::sqrt (pow (latvec.e31 , 2 ) + pow (latvec.e32 , 2 )
1217- + pow (latvec.e33 , 2 ));
1218- double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22
1219- + latvec.e13 * latvec.e23 );
1220- double cos12 = celldm12 / celldm1 / celldm2;
1221-
1222- double e21 = celldm2 * cos12;
1223- double e22 = celldm2 * std::sqrt (1.0 - cos12 * cos12);
1224-
1225- latvec.e11 = celldm1;
1226- latvec.e12 = 0.0 ;
1227- latvec.e13 = 0.0 ;
1228- latvec.e21 = e21 ;
1229- latvec.e22 = e22 ;
1230- latvec.e23 = 0.0 ;
1231- latvec.e31 = 0.0 ;
1232- latvec.e32 = 0.0 ;
1233- latvec.e33 = celldm3;
1234- } else if (latName == " bacm" ) // ibrav = 13
1235- {
1236- double celldm1 = std::abs (latvec.e11 );
1237- double celldm2 = std::sqrt (pow (latvec.e21 , 2 ) + pow (latvec.e22 , 2 )
1238- + pow (latvec.e23 , 2 ));
1239- double celldm3 = std::abs (latvec.e13 );
1240-
1241- double cos12 = latvec.e21 / celldm2;
1242- if (cos12 >= 1.0 ) {
1243- ModuleBase::WARNING_QUIT (" unitcell" , " wrong cos12!" );
1244- }
1245-
1246- double e21 = celldm2 * cos12;
1247- double e22 = celldm2 * std::sqrt (1.0 - cos12 * cos12);
1248-
1249- latvec.e11 = celldm1;
1250- latvec.e12 = 0.0 ;
1251- latvec.e13 = -celldm3;
1252- latvec.e21 = e21 ;
1253- latvec.e22 = e22 ;
1254- latvec.e23 = 0.0 ;
1255- latvec.e31 = celldm1;
1256- latvec.e32 = 0.0 ;
1257- latvec.e33 = celldm3;
1258- } else if (latName == " triclinic" ) // ibrav = 14
1259- {
1260- double celldm1 = std::sqrt (pow (latvec.e11 , 2 ) + pow (latvec.e12 , 2 )
1261- + pow (latvec.e13 , 2 ));
1262- double celldm2 = std::sqrt (pow (latvec.e21 , 2 ) + pow (latvec.e22 , 2 )
1263- + pow (latvec.e23 , 2 ));
1264- double celldm3 = std::sqrt (pow (latvec.e31 , 2 ) + pow (latvec.e32 , 2 )
1265- + pow (latvec.e33 , 2 ));
1266- double celldm12 = (latvec.e11 * latvec.e21 + latvec.e12 * latvec.e22
1267- + latvec.e13 * latvec.e23 );
1268- double cos12 = celldm12 / celldm1 / celldm2;
1269- double celldm13 = (latvec.e11 * latvec.e31 + latvec.e12 * latvec.e32
1270- + latvec.e13 * latvec.e33 );
1271- double cos13 = celldm13 / celldm1 / celldm3;
1272- double celldm23 = (latvec.e21 * latvec.e31 + latvec.e22 * latvec.e32
1273- + latvec.e23 * latvec.e33 );
1274- double cos23 = celldm23 / celldm2 / celldm3;
1275-
1276- double sin12 = std::sqrt (1.0 - cos12 * cos12);
1277- if (cos12 >= 1.0 ) {
1278- ModuleBase::WARNING_QUIT (" unitcell" , " wrong cos12!" );
1279- }
1280-
1281- latvec.e11 = celldm1;
1282- latvec.e12 = 0.0 ;
1283- latvec.e13 = 0.0 ;
1284- latvec.e21 = celldm2 * cos12;
1285- latvec.e22 = celldm2 * sin12;
1286- latvec.e23 = 0.0 ;
1287- latvec.e31 = celldm3 * cos13;
1288- latvec.e32 = celldm3 * (cos23 - cos13 * cos12) / sin12;
1289- double term = 1.0 + 2.0 * cos12 * cos13 * cos23 - cos12 * cos12
1290- - cos13 * cos13 - cos23 * cos23;
1291- term = sqrt (term) / sin12;
1292- latvec.e33 = celldm3 * term;
1293- } else {
1294- std::cout << " latname is : " << latName << std::endl;
1295- ModuleBase::WARNING_QUIT (" UnitCell::read_atom_species" ,
1296- " latname not supported!" );
1297- }
1298- }
12991023
13001024void UnitCell::compare_atom_labels (std::string label1, std::string label2) {
13011025 if (label1
0 commit comments