mirror of https://gitlab.com/QEF/q-e.git
259 lines
14 KiB
C
259 lines
14 KiB
C
/*
|
|
* Copyright (C) 2016 Quantum ESPRESSO group
|
|
* This file is distributed under the terms of the
|
|
* GNU General Public License. See the file `License'
|
|
* in the root directory of the present distribution,
|
|
* or http://www.gnu.org/copyleft/gpl.txt .
|
|
*/
|
|
|
|
/* ----------------------------------------------------------------------
|
|
This file contains auxiliary fuctions to interface QE with LAMMPS
|
|
|
|
Contributing authors: Axel Kohlmeyer (ICTP), Carlo Cavazzoni (CINECA)
|
|
------------------------------------------------------------------------- */
|
|
|
|
#include<stdio.h>
|
|
#include<stdlib.h>
|
|
#include<math.h>
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
/* Manage the atomic number */
|
|
|
|
#define QMMM_ISOTOPES 351
|
|
|
|
static const int FixQMMM_Z[QMMM_ISOTOPES] = {
|
|
1, 1, 1, 2, 2, 3, 3, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 8, 9,
|
|
10, 10, 10, 11, 12, 12, 12, 13, 14, 14, 14, 15, 16, 16, 16,
|
|
16, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 20, 20, 20,
|
|
21, 22, 22, 22, 22, 22, 23, 23, 24, 24, 24, 24, 25, 26, 26,
|
|
26, 26, 27, 28, 28, 28, 28, 28, 29, 29, 30, 30, 30, 30, 30,
|
|
31, 31, 32, 32, 32, 32, 32, 33, 34, 34, 34, 34, 34, 34, 35,
|
|
35, 36, 36, 36, 36, 36, 36, 37, 37, 38, 38, 38, 38, 39, 40,
|
|
40, 40, 40, 40, 41, 42, 42, 42, 42, 42, 42, 42, 43, 43, 43,
|
|
44, 44, 44, 44, 44, 44, 44, 45, 46, 46, 46, 46, 46, 46, 47,
|
|
47, 48, 48, 48, 48, 48, 48, 48, 48, 49, 49, 50, 50, 50, 50,
|
|
50, 50, 50, 50, 50, 50, 51, 51, 52, 52, 52, 52, 52, 52, 52,
|
|
52, 53, 54, 54, 54, 54, 54, 54, 54, 54, 54, 55, 56, 56, 56,
|
|
56, 56, 56, 56, 57, 57, 58, 58, 58, 58, 59, 60, 60, 60, 60,
|
|
60, 60, 60, 61, 61, 62, 62, 62, 62, 62, 62, 62, 63, 63, 64,
|
|
64, 64, 64, 64, 64, 64, 65, 66, 66, 66, 66, 66, 66, 66, 67,
|
|
68, 68, 68, 68, 68, 68, 69, 70, 70, 70, 70, 70, 70, 70, 71,
|
|
71, 72, 72, 72, 72, 72, 72, 73, 73, 74, 74, 74, 74, 74, 75,
|
|
75, 76, 76, 76, 76, 76, 76, 76, 77, 77, 78, 78, 78, 78, 78,
|
|
78, 79, 80, 80, 80, 80, 80, 80, 80, 81, 81, 82, 82, 82, 82,
|
|
83, 84, 84, 85, 85, 86, 86, 86, 87, 88, 88, 88, 88, 89, 90,
|
|
90, 91, 92, 92, 92, 92, 92, 93, 93, 94, 94, 94, 94, 94, 94,
|
|
95, 95, 96, 96, 96, 96, 96, 96, 97, 97, 98, 98, 98, 98, 99,
|
|
100, 101, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110,
|
|
111, 112, 113, 114, 115};
|
|
|
|
static const char FixQMMM_EL[QMMM_ISOTOPES][4] = {
|
|
"H", "H", "H", "He", "He", "Li", "Li", "Be", "B", "B",
|
|
"C", "C", "C", "N", "N", "O", "O", "O", "F", "Ne", "Ne",
|
|
"Ne", "Na", "Mg", "Mg", "Mg", "Al", "Si", "Si", "Si",
|
|
"P", "S", "S", "S", "S", "Cl", "Cl", "Ar", "Ar", "Ar",
|
|
"K", "K", "K", "Ca", "Ca", "Ca", "Ca", "Ca", "Ca", "Sc",
|
|
"Ti", "Ti", "Ti", "Ti", "Ti", "V", "V", "Cr", "Cr", "Cr",
|
|
"Cr", "Mn", "Fe", "Fe", "Fe", "Fe", "Co", "Ni", "Ni",
|
|
"Ni", "Ni", "Ni", "Cu", "Cu", "Zn", "Zn", "Zn", "Zn", "Zn",
|
|
"Ga", "Ga", "Ge", "Ge", "Ge", "Ge", "Ge", "As", "Se",
|
|
"Se", "Se", "Se", "Se", "Se", "Br", "Br", "Kr", "Kr",
|
|
"Kr", "Kr", "Kr", "Kr", "Rb", "Rb", "Sr", "Sr", "Sr",
|
|
"Sr", "Y", "Zr", "Zr", "Zr", "Zr", "Zr", "Nb", "Mo", "Mo",
|
|
"Mo", "Mo", "Mo", "Mo", "Mo", "Tc", "Tc", "Tc", "Ru", "Ru",
|
|
"Ru", "Ru", "Ru", "Ru", "Ru", "Rh", "Pd", "Pd", "Pd", "Pd",
|
|
"Pd", "Pd", "Ag", "Ag", "Cd", "Cd", "Cd", "Cd", "Cd", "Cd",
|
|
"Cd", "Cd", "In", "In", "Sn", "Sn", "Sn", "Sn", "Sn", "Sn",
|
|
"Sn", "Sn", "Sn", "Sn", "Sb", "Sb", "Te", "Te", "Te", "Te",
|
|
"Te", "Te", "Te", "Te", "I", "Xe", "Xe", "Xe", "Xe", "Xe",
|
|
"Xe", "Xe", "Xe", "Xe", "Cs", "Ba", "Ba", "Ba", "Ba", "Ba",
|
|
"Ba", "Ba", "La", "La", "Ce", "Ce", "Ce", "Ce", "Pr", "Nd",
|
|
"Nd", "Nd", "Nd", "Nd", "Nd", "Nd", "Pm", "Pm", "Sm", "Sm",
|
|
"Sm", "Sm", "Sm", "Sm", "Sm", "Eu", "Eu", "Gd", "Gd", "Gd",
|
|
"Gd", "Gd", "Gd", "Gd", "Tb", "Dy", "Dy", "Dy", "Dy", "Dy",
|
|
"Dy", "Dy", "Ho", "Er", "Er", "Er", "Er", "Er", "Er", "Tm",
|
|
"Yb", "Yb", "Yb", "Yb", "Yb", "Yb", "Yb", "Lu", "Lu", "Hf",
|
|
"Hf", "Hf", "Hf", "Hf", "Hf", "Ta", "Ta", "W", "W", "W", "W",
|
|
"W", "Re", "Re", "Os", "Os", "Os", "Os", "Os", "Os", "Os",
|
|
"Ir", "Ir", "Pt", "Pt", "Pt", "Pt", "Pt", "Pt", "Au", "Hg",
|
|
"Hg", "Hg", "Hg", "Hg", "Hg", "Hg", "Tl", "Tl", "Pb", "Pb",
|
|
"Pb", "Pb", "Bi", "Po", "Po", "At", "At", "Rn", "Rn", "Rn",
|
|
"Fr", "Ra", "Ra", "Ra", "Ra", "Ac", "Th", "Th", "Pa", "U",
|
|
"U", "U", "U", "U", "Np", "Np", "Pu", "Pu", "Pu", "Pu", "Pu",
|
|
"Pu", "Am", "Am", "Cm", "Cm", "Cm", "Cm", "Cm", "Cm", "Bk",
|
|
"Bk", "Cf", "Cf", "Cf", "Cf", "Es", "Fm", "Md", "Md", "No",
|
|
"Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", "Cn",
|
|
"Uut", "Uuq", "Uup"};
|
|
|
|
static const int FixQMMM_A[QMMM_ISOTOPES] = {
|
|
1, 2, 3, 3, 4, 6, 7, 9, 10, 11, 12, 13, 14, 14, 15, 16, 17, 18,
|
|
19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
|
|
36, 35, 37, 36, 38, 40, 39, 40, 41, 40, 42, 43, 44, 46, 48, 45,
|
|
46, 47, 48, 49, 50, 50, 51, 50, 52, 53, 54, 55, 54, 56, 57, 58,
|
|
59, 58, 60, 61, 62, 64, 63, 65, 64, 66, 67, 68, 70, 69, 71, 70,
|
|
72, 73, 74, 76, 75, 74, 76, 77, 78, 80, 82, 79, 81, 78, 80, 82,
|
|
83, 84, 86, 85, 87, 84, 86, 87, 88, 89, 90, 91, 92, 94, 96, 93,
|
|
92, 94, 95, 96, 97, 98, 100, 97, 98, 99, 96, 98, 99, 100, 101,
|
|
102, 104, 103, 102, 104, 105, 106, 108, 110, 107, 109, 106, 108,
|
|
110, 111, 112, 113, 114, 116, 113, 115, 112, 114, 115, 116, 117,
|
|
118, 119, 120, 122, 124, 121, 123, 120, 122, 123, 124, 125, 126,
|
|
128, 130, 127, 124, 126, 128, 129, 130, 131, 132, 134, 136, 133,
|
|
130, 132, 134, 135, 136, 137, 138, 138, 139, 136, 138, 140, 142,
|
|
141, 142, 143, 144, 145, 146, 148, 150, 145, 147, 144, 147, 148,
|
|
149, 150, 152, 154, 151, 153, 152, 154, 155, 156, 157, 158, 160,
|
|
159, 156, 158, 160, 161, 162, 163, 164, 165, 162, 164, 166, 167,
|
|
168, 170, 169, 168, 170, 171, 172, 173, 174, 176, 175, 176, 174,
|
|
176, 177, 178, 179, 180, 180, 181, 180, 182, 183, 184, 186, 185,
|
|
187, 184, 186, 187, 188, 189, 190, 192, 191, 193, 190, 192, 194,
|
|
195, 196, 198, 197, 196, 198, 199, 200, 201, 202, 204, 203, 205,
|
|
204, 206, 207, 208, 209, 209, 210, 210, 211, 211, 220, 222, 223,
|
|
223, 224, 226, 228, 227, 230, 232, 231, 233, 234, 235, 236, 238,
|
|
236, 237, 238, 239, 240, 241, 242, 244, 241, 243, 243, 244, 245,
|
|
246, 247, 248, 247, 249, 249, 250, 251, 252, 252, 257, 258, 260,
|
|
259, 262, 265, 268, 271, 272, 270, 276, 281, 280, 285, 284, 289,
|
|
288};
|
|
|
|
static const double FixQMMM_MASS[QMMM_ISOTOPES] = {
|
|
1.00782503207, 2.0141017778, 3.0160492777, 3.0160293191, 4.00260325415, 6.015122795, 7.01600455, 9.0121822,
|
|
10.0129370, 11.0093054, 12.0000000, 13.0033548378, 14.003241989, 14.0030740048, 15.0001088982, 15.99491461956,
|
|
16.99913170, 17.9991610, 18.99840322, 19.9924401754, 20.99384668, 21.991385114, 22.9897692809, 23.985041700,
|
|
24.98583692, 25.982592929, 26.98153863, 27.9769265325, 28.976494700, 29.97377017, 30.97376163, 31.97207100,
|
|
32.97145876, 33.96786690, 35.96708076, 34.96885268, 36.96590259, 35.967545106, 37.9627324, 39.9623831225,
|
|
38.96370668, 39.96399848, 40.96182576, 39.96259098, 41.95861801, 42.9587666, 43.9554818, 45.9536926,
|
|
47.952534, 44.9559119, 45.9526316, 46.9517631, 47.9479463, 48.9478700, 49.9447912, 49.9471585,
|
|
50.9439595, 49.9460442, 51.9405075, 52.9406494, 53.9388804, 54.9380451, 53.9396105, 55.9349375,
|
|
56.9353940, 57.9332756, 58.9331950, 57.9353429, 59.9307864, 60.9310560, 61.9283451, 63.9279660,
|
|
62.9295975, 64.9277895, 63.9291422, 65.9260334, 66.9271273, 67.9248442, 69.9253193, 68.9255736,
|
|
70.9247013, 69.9242474, 71.9220758, 72.9234589, 73.9211778, 75.9214026, 74.9215965, 73.9224764,
|
|
75.9192136, 76.9199140, 77.9173091, 79.9165213, 81.9166994, 78.9183371, 80.9162906, 77.9203648,
|
|
79.9163790, 81.9134836, 82.914136, 83.911507, 85.91061073, 84.911789738, 86.909180527, 83.913425,
|
|
85.9092602, 86.9088771, 87.9056121, 88.9058483, 89.9047044, 90.9056458, 91.9050408, 93.9063152,
|
|
95.9082734, 92.9063781, 91.906811, 93.9050883, 94.9058421, 95.9046795, 96.9060215, 97.9054082,
|
|
99.907477, 96.906365, 97.907216, 98.9062547, 95.907598, 97.905287, 98.9059393, 99.9042195,
|
|
100.9055821, 101.9043493, 103.905433, 102.905504, 101.905609, 103.904036, 104.905085, 105.903486,
|
|
107.903892, 109.905153, 106.905097, 108.904752, 105.906459, 107.904184, 109.9030021, 110.9041781,
|
|
111.9027578, 112.9044017, 113.9033585, 115.904756, 112.904058, 114.903878, 111.904818, 113.902779,
|
|
114.903342, 115.901741, 116.902952, 117.901603, 118.903308, 119.9021947, 121.9034390, 123.9052739,
|
|
120.9038157, 122.9042140, 119.904020, 121.9030439, 122.9042700, 123.9028179, 124.9044307, 125.9033117,
|
|
127.9044631, 129.9062244, 126.904473, 123.9058930, 125.904274, 127.9035313, 128.9047794, 129.9035080,
|
|
130.9050824, 131.9041535, 133.9053945, 135.907219, 132.905451933, 129.9063208, 131.9050613, 133.9045084,
|
|
134.9056886, 135.9045759, 136.9058274, 137.9052472, 137.907112, 138.9063533, 135.907172, 137.905991,
|
|
139.9054387, 141.909244, 140.9076528, 141.9077233, 142.9098143, 143.9100873, 144.9125736, 145.9131169,
|
|
147.916893, 149.920891, 144.912749, 146.9151385, 143.911999, 146.9148979, 147.9148227, 148.9171847,
|
|
149.9172755, 151.9197324, 153.9222093, 150.9198502, 152.9212303, 151.9197910, 153.9208656, 154.9226220,
|
|
155.9221227, 156.9239601, 157.9241039, 159.9270541, 158.9253468, 155.924283, 157.924409, 159.9251975,
|
|
160.9269334, 161.9267984, 162.9287312, 163.9291748, 164.9303221, 161.928778, 163.929200, 165.9302931,
|
|
166.9320482, 167.9323702, 169.9354643, 168.9342133, 167.933897, 169.9347618, 170.9363258, 171.9363815,
|
|
172.9382108, 173.9388621, 175.9425717, 174.9407718, 175.9426863, 173.940046, 175.9414086, 176.9432207,
|
|
177.9436988, 178.9458161, 179.9465500, 179.9474648, 180.9479958, 179.946704, 181.9482042, 182.9502230,
|
|
183.9509312, 185.9543641, 184.9529550, 186.9557531, 183.9524891, 185.9538382, 186.9557505, 187.9558382,
|
|
188.9581475, 189.9584470, 191.9614807, 190.9605940, 192.9629264, 189.959932, 191.9610380, 193.9626803,
|
|
194.9647911, 195.9649515, 197.967893, 196.9665687, 195.965833, 197.9667690, 198.9682799, 199.9683260,
|
|
200.9703023, 201.9706430, 203.9734939, 202.9723442, 204.9744275, 203.9730436, 205.9744653, 206.9758969,
|
|
207.9766521, 208.9803987, 208.9824304, 209.9828737, 209.987148, 210.9874963, 210.990601, 220.0113940,
|
|
222.0175777, 223.0197359, 223.0185022, 224.0202118, 226.0254098, 228.0310703, 227.0277521, 230.0331338,
|
|
232.0380553, 231.0358840, 233.0396352, 234.0409521, 235.0439299, 236.0455680, 238.0507882, 236.046570,
|
|
237.0481734, 238.0495599, 239.0521634, 240.0538135, 241.0568515, 242.0587426, 244.064204, 241.0568291,
|
|
243.0613811, 243.0613891, 244.0627526, 245.0654912, 246.0672237, 247.070354, 248.072349, 247.070307,
|
|
249.0749867, 249.0748535, 250.0764061, 251.079587, 252.081626, 252.082980, 257.095105, 258.098431,
|
|
260.10365, 259.10103, 262.10963, 265.11670, 268.12545, 271.13347, 272.13803, 270.13465,
|
|
276.15116, 281.16206, 280.16447, 285.17411, 284.17808, 289.18728, 288.19249};
|
|
|
|
/*
|
|
* This function matches the element which has the least absolute
|
|
* difference between the masses. delta is set to the difference
|
|
* between the mass provided and the best match.
|
|
*
|
|
* THIS FUNCTION RELIES ON THE CORRECT ORDER OF THE DATA TABLES
|
|
* IN ORDER TO SKIP ISOTOPES FROM THE MATCH
|
|
*
|
|
* */
|
|
int match_element (double mass, int search_isotopes, double * delta)
|
|
{
|
|
int i;
|
|
int bestcandidate;
|
|
double diff, mindiff;
|
|
int lastz;
|
|
|
|
mindiff = 1e6;
|
|
bestcandidate = -1;
|
|
lastz = -1;
|
|
for(i=0; i<QMMM_ISOTOPES; i++) {
|
|
if(!search_isotopes && lastz == FixQMMM_Z[i])
|
|
continue;
|
|
diff = fabs(FixQMMM_MASS[i] - mass);
|
|
if(diff < mindiff) {
|
|
mindiff = diff;
|
|
bestcandidate = i;
|
|
}
|
|
lastz = FixQMMM_Z[i];
|
|
}
|
|
*delta = mindiff;
|
|
return bestcandidate;
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
//! Atomic radii of all atoms. Master/QM only
|
|
//
|
|
//
|
|
|
|
/*
|
|
* Source: wikipedia
|
|
* [ http://en.wikipedia.org/wiki/Atomic_radii_of_the_elements_(data_page) ]
|
|
*
|
|
* Atomic radii for the various elements in picometers.
|
|
*
|
|
* A value of -1 has the meaning of "data not available".
|
|
*
|
|
* */
|
|
|
|
/// Number of elements in the arrays
|
|
#define EC_ELEMENTS 116
|
|
|
|
static const double ec_r_covalent[EC_ELEMENTS+1] = {
|
|
-1 /* Z=0 */,
|
|
0.38, 0.32, 1.34, 0.90, 0.82, 0.77, 0.75, 0.73, 0.71, 0.69, 1.54, 1.30, 1.18, 1.11, 1.06, 1.02, 0.99, 0.97, 1.96, 1.74,
|
|
1.44, 1.36, 1.25, 1.27, 1.39, 1.25, 1.26, 1.21, 1.38, 1.31, 1.26, 1.22, 1.19, 1.16, 1.14, 1.10, 2.11, 1.92, 1.62, 1.48,
|
|
1.37, 1.45, 1.56, 1.26, 1.35, 1.31, 1.53, 1.48, 1.44, 1.41, 1.38, 1.35, 1.33, 1.30, 2.25, 1.98, 1.69, -1, -1, -1,
|
|
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1.60, 1.50, 1.38, 1.46, 1.59, 1.28, 1.37, 1.28, 1.44, 1.49,
|
|
1.48, 1.47, 1.46, -1, -1, 1.45, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
|
|
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1
|
|
};
|
|
|
|
int ec_fill_radii
|
|
( double *aradii, int *ndata, double *mass, int *type, int *ntypes, int *screen )
|
|
{
|
|
int i, el;
|
|
int *type2element;
|
|
|
|
type2element = (int *) malloc( sizeof(int) * *ntypes+1 );
|
|
type2element[0] = 0;
|
|
|
|
for (i=1; i <= *ntypes; ++i) {
|
|
double delta;
|
|
el = match_element(mass[i], 0, &delta);
|
|
if(*screen) fprintf(stdout,"FixQMMM: type %2d (mass: %8g) matches %2s with:"
|
|
" Z = %-3d A = %-3d r_cov = %5.2f (error = %-8.2g -> %-.2g%%)\n",
|
|
i, mass[i], FixQMMM_EL[el], FixQMMM_Z[el], FixQMMM_A[el],
|
|
ec_r_covalent[FixQMMM_Z[el]], delta, delta/mass[i] * 100.0);
|
|
type2element[i] = FixQMMM_Z[el];
|
|
}
|
|
|
|
for (i=0; i<*ndata; i++) {
|
|
el = type2element[type[i]];
|
|
if (el < 0 || el > EC_ELEMENTS) {
|
|
fprintf(stderr,"Unable to find element Z=%d in table of covalent radii", el);
|
|
exit(1);
|
|
}
|
|
aradii[i] = ec_r_covalent[el];
|
|
if (ec_r_covalent[el] < 0.0) {
|
|
fprintf(stderr, "Covalent radius for atom of element Z=%d not availabe", el);
|
|
exit(1);
|
|
}
|
|
}
|
|
|
|
free( (void *) type2element );
|
|
return 0;
|
|
}
|
|
|