Line data    Source code

       1             : /*
       2             :  * Copyright (c) 2019 Sebastian Weber, Henri Menke, Johannes Block. All rights reserved.
       3             :  *
       4             :  * This file is part of the pairinteraction library.
       5             :  *
       6             :  * The pairinteraction library is free software: you can redistribute it and/or modify
       7             :  * it under the terms of the GNU Lesser General Public License as published by
       8             :  * the Free Software Foundation, either version 3 of the License, or
       9             :  * (at your option) any later version.
      10             :  *
      11             :  * The pairinteraction library is distributed in the hope that it will be useful,
      12             :  * but WITHOUT ANY WARRANTY; without even the implied warranty of
      13             :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      14             :  * GNU Lesser General Public License for more details.
      15             :  *
      16             :  * You should have received a copy of the GNU Lesser General Public License
      17             :  * along with the pairinteraction library. If not, see <http://www.gnu.org/licenses/>.
      18             :  */
      19             : 
      20             : #include "GreenTensor.hpp"
      21             : 
      22             : #include "EigenCompat.hpp"
      23             : #include <Eigen/Core>
      24             : #include <unsupported/Eigen/CXX11/Tensor>
      25             : 
      26             : #include <limits>
      27             : 
      28             : using TensorType = GreenTensor::TensorType;
      29             : 
      30          42 : GreenTensor::GreenTensor(double x, double y, double z)
      31          42 :     : x(x), y(y), z(z), zA(std::numeric_limits<double>::max()),
      32          42 :       zB(std::numeric_limits<double>::max()), dd_tensor_calculated(false),
      33          42 :       dq_tensor_calculated(false), qd_tensor_calculated(false) {}
      34             : 
      35          37 : void GreenTensor::addSurface(double d) {
      36          37 :     if (y != 0) {
      37           0 :         throw std::runtime_error("The atoms must be in the xz-plane if a surface is present");
      38             :     }
      39             : 
      40          37 :     double angle = std::atan(x / z);
      41             : 
      42          37 :     zA = d - z * std::sin(angle) / 2.;
      43          37 :     zB = d + z * std::sin(angle) / 2.;
      44             : 
      45          37 :     if (zA < 0 || zB < 0) {
      46           0 :         throw std::runtime_error(
      47           0 :             "zA or zB < 0. One of the atoms is inside the plate. Plate is half-space z < 0.");
      48             :     }
      49             : 
      50          37 :     dd_tensor_calculated = false;
      51          37 :     dq_tensor_calculated = false;
      52          37 :     qd_tensor_calculated = false;
      53          37 : }
      54             : 
      55          42 : const Eigen::Matrix3<double> &GreenTensor::getDDTensor() {
      56          42 :     if (!dd_tensor_calculated) {
      57          42 :         dd_tensor = getDDTensorVacuum(x, y, z);
      58          42 :         if (zA != std::numeric_limits<double>::max()) {
      59          37 :             dd_tensor += getDDTensorPlate(x, zA, zB);
      60             :         }
      61          42 :         dd_tensor_calculated = true;
      62             :     }
      63             : 
      64          42 :     return dd_tensor;
      65             : }
      66             : 
      67           1 : const TensorType &GreenTensor::getDQTensor() {
      68           1 :     if (!dq_tensor_calculated) {
      69           1 :         dq_tensor = getDQTensorVacuum(x, y, z);
      70           1 :         if (zA != std::numeric_limits<double>::max()) {
      71           0 :             dq_tensor += getDQTensorPlate(x, zA, zB);
      72             :         }
      73           1 :         dq_tensor_calculated = true;
      74             :     }
      75           1 :     return dq_tensor;
      76             : }
      77             : 
      78           1 : const TensorType &GreenTensor::getQDTensor() {
      79           1 :     if (!qd_tensor_calculated) {
      80           1 :         qd_tensor = getQDTensorVacuum(x, y, z);
      81           1 :         if (zA != std::numeric_limits<double>::max()) {
      82           0 :             qd_tensor += getQDTensorPlate(x, zA, zB);
      83             :         }
      84           1 :         qd_tensor_calculated = true;
      85             :     }
      86           1 :     return qd_tensor;
      87             : }
      88             : 
      89          42 : Eigen::Matrix3<double> GreenTensor::getDDTensorVacuum(double x, double y, double z) const {
      90             :     // Build distance vector
      91          42 :     Eigen::Matrix<double, 3, 1> distance;
      92          42 :     distance << x, y, z;
      93             : 
      94             :     // Construct Green tensor
      95             :     Eigen::Matrix3<double> vacuum_tensor =
      96          42 :         -Eigen::MatrixXd::Identity(3, 3) / std::pow(distance.norm(), 3.) +
      97          84 :         3. * distance * distance.transpose() / std::pow(distance.norm(), 5.);
      98             : 
      99          84 :     return vacuum_tensor;
     100             : }
     101             : 
     102          37 : Eigen::Matrix3<double> GreenTensor::getDDTensorPlate(double x, double zA, double zB) const {
     103             :     // Calculate distances to mirror dipole
     104          37 :     double zp = zA + zB;
     105          37 :     double rp = std::sqrt(x * x + zp * zp);
     106             : 
     107             :     // Construct Green tensor
     108          37 :     Eigen::Matrix3<double> plate_tensor_second_matrix;
     109          37 :     plate_tensor_second_matrix << x * x, 0., -x * zp, 0., 0., 0., x * zp, 0., x * x;
     110             : 
     111             :     Eigen::Matrix3<double> plate_tensor =
     112          37 :         Eigen::Vector3<double>({1., 1., 2.}).asDiagonal().toDenseMatrix() / std::pow(rp, 3.) -
     113          74 :         3. * plate_tensor_second_matrix / std::pow(rp, 5);
     114             : 
     115          74 :     return plate_tensor;
     116             : }
     117             : 
     118           1 : TensorType GreenTensor::getDQTensorVacuum(double x, double y, double z) const {
     119           1 :     TensorType vacuum_tensor;
     120           1 :     vacuum_tensor.setZero();
     121             : 
     122           1 :     Eigen::Matrix3<double> Eye = Eigen::Matrix3<double>::Identity(3, 3);
     123           1 :     Eigen::Matrix<double, 3, 1> distance;
     124           1 :     distance << x, y, z;
     125           1 :     double dist = std::sqrt(x * x + y * y + z * z);
     126             : 
     127           4 :     for (int i = 0; i < 3; i++) {
     128          12 :         for (int j = 0; j < 3; j++) {
     129          36 :             for (int k = 0; k < 3; k++) {
     130          27 :                 vacuum_tensor(i, j, k) += -3 / std::pow(dist, 5.) *
     131          27 :                         (Eye(i, j) * distance(k) + distance(i) * Eye(j, k) +
     132          27 :                          distance(j) * Eye(i, k)) +
     133          27 :                     15. / std::pow(dist, 7.) * distance(i) * distance(j) * distance(k);
     134             :             }
     135             :         }
     136             :     }
     137             : 
     138           2 :     return vacuum_tensor;
     139             : }
     140             : 
     141           1 : TensorType GreenTensor::getQDTensorVacuum(double x, double y, double z) const {
     142           1 :     TensorType vacuum_tensor;
     143           1 :     vacuum_tensor.setZero();
     144             : 
     145           1 :     Eigen::Matrix3<double> Eye = Eigen::Matrix3<double>::Identity(3, 3);
     146           1 :     Eigen::Matrix<double, 3, 1> distance;
     147           1 :     distance << x, y, z;
     148           1 :     double dist = std::sqrt(x * x + y * y + z * z);
     149             : 
     150           4 :     for (int i = 0; i < 3; i++) {
     151          12 :         for (int j = 0; j < 3; j++) {
     152          36 :             for (int k = 0; k < 3; k++) {
     153             : 
     154          27 :                 vacuum_tensor(i, j, k) += 3 / std::pow(dist, 5.) *
     155          27 :                         (distance(i) * Eye(j, k) + Eye(i, j) * distance(k) +
     156          27 :                          Eye(i, k) * distance(j)) -
     157          27 :                     15. / std::pow(dist, 7.) * distance(i) * distance(j) * distance(k);
     158             :             }
     159             :         }
     160             :     }
     161             : 
     162           2 :     return vacuum_tensor;
     163             : }
     164             : 
     165           0 : TensorType GreenTensor::getDQTensorPlate(double x, double zA, double zB) const {
     166           0 :     TensorType plate_tensor;
     167           0 :     plate_tensor.setZero();
     168             : 
     169           0 :     double zp = zA + zB;
     170           0 :     Eigen::Matrix<double, 3, 1> distanceplus;
     171           0 :     distanceplus << -x, 0., zp;
     172           0 :     double rp = distanceplus.norm();
     173             : 
     174             :     Eigen::Matrix3<double> plate_tensor_first_matrix =
     175           0 :         Eigen::Vector3<double>({1., 1., 2.}).asDiagonal().toDenseMatrix();
     176           0 :     Eigen::Matrix3<double> plate_tensor_second_matrix;
     177           0 :     plate_tensor_second_matrix << x * x, 0., -x * zp, 0., 0., 0., x * zp, 0., x * x;
     178             : 
     179           0 :     TensorType plate_tensor_second_gradient;
     180           0 :     plate_tensor_second_gradient.setZero();
     181           0 :     plate_tensor_second_gradient(0, 0, 0) = -2. * x;
     182           0 :     plate_tensor_second_gradient(2, 2, 0) = -2. * x;
     183           0 :     plate_tensor_second_gradient(0, 2, 0) = zp;
     184           0 :     plate_tensor_second_gradient(2, 0, 0) = -zp;
     185           0 :     plate_tensor_second_gradient(0, 2, 2) = -x;
     186           0 :     plate_tensor_second_gradient(2, 0, 2) = x;
     187             : 
     188           0 :     for (int i = 0; i < 3; i++) {
     189           0 :         for (int j = 0; j < 3; j++) {
     190           0 :             for (int k = 0; k < 3; k++) {
     191           0 :                 plate_tensor(i, j, k) +=
     192           0 :                     -3. * plate_tensor_first_matrix(i, j) * distanceplus(k) / std::pow(rp, 5.);
     193           0 :                 plate_tensor(i, j, k) +=
     194           0 :                     15. * plate_tensor_second_matrix(i, j) * distanceplus(k) / std::pow(rp, 7.);
     195           0 :                 plate_tensor(i, j, k) +=
     196           0 :                     -3. * plate_tensor_second_gradient(i, j, k) / std::pow(rp, 5.);
     197             :             }
     198             :         }
     199             :     }
     200             : 
     201           0 :     return plate_tensor;
     202             : }
     203             : 
     204           0 : TensorType GreenTensor::getQDTensorPlate(double x, double zA, double zB) const {
     205           0 :     TensorType plate_tensor;
     206           0 :     plate_tensor.setZero();
     207             : 
     208           0 :     double zp = zA + zB;
     209           0 :     Eigen::Matrix<double, 3, 1> distanceplus;
     210           0 :     distanceplus << x, 0., zp;
     211           0 :     double rp = distanceplus.norm();
     212             : 
     213             :     Eigen::Matrix3<double> plate_tensor_first_matrix =
     214           0 :         Eigen::Vector3<double>({1., 1., 2.}).asDiagonal().toDenseMatrix();
     215           0 :     Eigen::Matrix3<double> plate_tensor_second_matrix;
     216           0 :     plate_tensor_second_matrix << x * x, 0., -x * zp, 0., 0., 0., x * zp, 0., x * x;
     217             : 
     218           0 :     TensorType gradient_plate_tensor_second;
     219           0 :     gradient_plate_tensor_second.setZero();
     220           0 :     gradient_plate_tensor_second(0, 0, 0) = 2. * x;
     221           0 :     gradient_plate_tensor_second(0, 2, 2) = 2. * x;
     222           0 :     gradient_plate_tensor_second(0, 0, 2) = -zp;
     223           0 :     gradient_plate_tensor_second(0, 2, 0) = zp;
     224           0 :     gradient_plate_tensor_second(2, 0, 2) = -x;
     225           0 :     gradient_plate_tensor_second(2, 2, 0) = x;
     226             : 
     227           0 :     for (int i = 0; i < 3; i++) {
     228           0 :         for (int j = 0; j < 3; j++) {
     229           0 :             for (int k = 0; k < 3; k++) {
     230           0 :                 plate_tensor(i, j, k) +=
     231           0 :                     -3. * distanceplus(i) * plate_tensor_first_matrix(j, k) / std::pow(rp, 5.);
     232           0 :                 plate_tensor(i, j, k) +=
     233           0 :                     15. * distanceplus(i) * plate_tensor_second_matrix(j, k) / std::pow(rp, 7.);
     234           0 :                 plate_tensor(i, j, k) +=
     235           0 :                     -3. * gradient_plate_tensor_second(i, j, k) / std::pow(rp, 5.);
     236             :             }
     237             :         }
     238             :     }
     239             : 
     240           0 :     return plate_tensor;
     241             : }