diff --git a/src/amoeba.F b/src/amoeba.F deleted file mode 100644 index cc56243..0000000 --- a/src/amoeba.F +++ /dev/null @@ -1,120 +0,0 @@ -!-----------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart ! -!-----------------------------------------------------------------------------! - -MODULE amoeba - - USE kinds, ONLY : dbl - USE nrutil, ONLY : assert_eq, imaxloc, iminloc, nrerror, swap - - IMPLICIT NONE - - PRIVATE - PUBLIC :: amoeba_evaluate - -CONTAINS - -!-----------------------------------------------------------------------------! -! NUMERICAL RECIPIES SUBROUTINE AMOEBA ADAPTED FOR USE IN CP2K ! -!-----------------------------------------------------------------------------! - -SUBROUTINE amoeba_evaluate ( p, y, ftol, rtol, func, itmax ) - - IMPLICIT NONE - -! Arguments - REAL (dbl), DIMENSION (:,:), INTENT (INOUT) :: p - REAL (dbl), DIMENSION (:), INTENT (INOUT) :: y - REAL (dbl), INTENT (IN) :: ftol - REAL (dbl), INTENT ( OUT ) :: rtol - INTEGER, INTENT ( IN ) :: itmax - -! Interface - INTERFACE - FUNCTION func ( x ) - USE kinds, ONLY : dbl - IMPLICIT NONE - REAL ( dbl ), DIMENSION ( : ), INTENT ( IN ) :: x - REAL ( dbl ) :: func - END FUNCTION func - END INTERFACE - - INTEGER :: ihi, ndim - REAL (dbl), DIMENSION (size(p,2)) :: psum - - CALL amoeba_private() - -!------------------------------------------------------------------------------ - -CONTAINS - - SUBROUTINE amoeba_private - IMPLICIT NONE - INTEGER :: i, ilo, inhi - INTEGER :: iter - REAL (dbl) :: ysave, ytry, ytmp - - ndim = assert_eq(size(p,2),size(p,1)-1,size(y)-1,'amoeba') - iter = 0 - psum(:) = sum(p(:,:),dim=1) - DO - ilo = iminloc(y(:)) - ihi = imaxloc(y(:)) - ytmp = y(ihi) - y(ihi) = y(ilo) - inhi = imaxloc(y(:)) - y(ihi) = ytmp - rtol = 2.0_dbl*abs(y(ihi)-y(ilo))/(abs(y(ihi))+abs(y(ilo))) - IF (rtol=itmax) RETURN - ytry = amotry(-1.0_dbl) - iter = iter + 1 - IF (ytry<=y(ilo)) THEN - ytry = amotry(2.0_dbl) - iter = iter + 1 - ELSE IF (ytry>=y(inhi)) THEN - ysave = y(ihi) - ytry = amotry(0.5_dbl) - iter = iter + 1 - IF (ytry>=ysave) THEN - p(:,:) = 0.5_dbl*(p(:,:)+spread(p(ilo,:),1,size(p,1))) - DO i = 1, ndim + 1 - IF (i/=ilo) y(i) = func(p(i,:)) - END DO - iter = iter + ndim - psum(:) = sum(p(:,:),dim=1) - END IF - END IF - END DO - END SUBROUTINE amoeba_private - -!BL - FUNCTION amotry(fac) - IMPLICIT NONE - REAL (dbl), INTENT (IN) :: fac - REAL (dbl) :: amotry - REAL (dbl) :: fac1, fac2, ytry - REAL (dbl), DIMENSION (size(p,2)) :: ptry - - fac1 = (1.0_dbl-fac)/ndim - fac2 = fac1 - fac - ptry(:) = psum(:)*fac1 - p(ihi,:)*fac2 - ytry = func(ptry) - IF (ytry +#include + +#include "amoeba.h" + +namespace amoeba { + +/*---------------------------------------------------------------------------*/ +/* NUMERICAL RECIPIES SUBROUTINE AMOEBA ADAPTED FOR USE IN CP2K */ +/*---------------------------------------------------------------------------*/ + + template + void amoeba_evaluate(std::vector>& p, std::vector& y, + double ftol, double& rtol, Func func, int itmax) { + int ndim = p[0].size(); + std::vector psum(ndim); + + amoeba_private(); + + //----------------------------------------------------------------------- + + void amoeba_private() { + int i, ilo, inhi; + int iter = 0; + double ysave, ytry, ytmp; + + ndim = p[0].size(); + psum.resize(ndim); + std::fill(psum.begin(), psum.end(), 0.0); + + do { + ilo = std::min_element(y.begin(), y.end()) - y.begin(); + ihi = std::max_element(y.begin(), y.end()) - y.begin(); + ytmp = y[ihi]; + y[ihi] = y[ilo]; + inhi = std::max_element(y.begin(), y.end()) - y.begin(); + y[ihi] = ytmp; + rtol = 2.0 * std::abs(y[ihi] - y[ilo]) / (std::abs(y[ihi]) + std::abs(y[ilo])); + if (rtol < ftol) { + std::swap(y[0], y[ilo]); + std::swap(p[0], p[ilo]); + return; + } + if (iter >= itmax) + return; + ytry = amotry(-1.0); + iter++; + if (ytry <= y[ilo]) { + ytry = amotry(2.0); + iter++; + } else if (ytry >= y[inhi]) { + ysave = y[ihi]; + ytry = amotry(0.5); + iter++; + if (ytry >= ysave) { + for (int i = 0; i < ndim + 1; i++) { + if (i != ilo) + y[i] = func(p[i]); + } + iter += ndim; + for (int j = 0; j < ndim; j++) { + psum[j] = 0.0; + for (int i = 0; i < ndim + 1; i++) { + psum[j] += p[i][j]; + } + } + } + } + } while (true); + } + + double amotry(double fac) { + double fac1 = (1.0 - fac) / ndim; + double fac2 = fac1 - fac; + std::vector ptry(ndim); + for (int j = 0; j < ndim; j++) { + ptry[j] = psum[j] * fac1 - p[ihi][j] * fac2; + } + double ytry = func(ptry); + if (ytry < y[ihi]) { + y[ihi] = ytry; + for (int j = 0; j < ndim; j++) { + psum[j] = psum[j] - p[ihi][j] + ptry[j]; + p[ihi][j] = ptry[j]; + } + } + return ytry; + } + } + +} // namespace amoeba diff --git a/src/amoeba.h b/src/amoeba.h new file mode 100644 index 0000000..9a33cfc --- /dev/null +++ b/src/amoeba.h @@ -0,0 +1,10 @@ +#ifndef _AMOEBA_H +#define _AMOEBA_H + +#include + +template +void amoeba_evaluate(std::vector>&, std::vector&, + double, double&, Func, int); + +#endif \ No newline at end of file diff --git a/src/atoms_input.F b/src/atoms_input.F deleted file mode 100644 index c04e912..0000000 --- a/src/atoms_input.F +++ /dev/null @@ -1,246 +0,0 @@ -!-----------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright (C) 2000 CP2K developers group ! -!-----------------------------------------------------------------------------! - -MODULE atoms_input - - USE global_types, ONLY : global_environment_type - USE kinds, ONLY : dbl - USE mp, ONLY : mp_bcast - USE parser, ONLY : parser_init, parser_end, read_line, test_next, & - cfield, p_error, get_real, get_int - USE stop_program, ONLY : stop_prg, stop_memory - USE string_utilities, ONLY : xstring, uppercase - USE util, ONLY : get_unit - - IMPLICIT NONE - - PRIVATE - PUBLIC :: read_coord_vel, system_type - - TYPE system_type - REAL ( dbl ) :: box(3,3) - CHARACTER ( LEN = 3 ) :: ptype - CHARACTER ( LEN = 20 ) :: rtype - INTEGER :: n - REAL ( dbl ), POINTER :: c ( :, : ) - REAL ( dbl ), POINTER :: v ( :, : ) - END TYPE system_type - -CONTAINS - -!!>---------------------------------------------------------------------------! -!! ! -!! OPTIONS: either read the section &atoms in the input file or ! -!! read the coordinates from the file project_name.dat ! -!! SECTION: &atoms ... &end ! -!! file: file_name ! -!! cell: b11 b12 b13 & ! -!! b21 b22 b23 & ! -!! b31 b32 b33 ! -!! ! -!! ! -!! ! -!!<---------------------------------------------------------------------------! - -SUBROUTINE read_coord_vel ( atype, filen, globenv ) - - IMPLICIT NONE - -! Arguments - TYPE ( system_type ), INTENT ( INOUT ) :: atype - CHARACTER ( LEN = * ), INTENT ( INOUT ) :: filen - TYPE ( global_environment_type ), INTENT ( IN ) :: globenv - -! Locals - INTEGER :: ierror, ilen, iw, source, allgrp, ia, ie, i, ios - CHARACTER ( LEN = 6 ) :: string - CHARACTER ( LEN = 5 ) :: label - -!------------------------------------------------------------------------------ - - iw = globenv % scr - -!..defaults -!..parse the input section - label = '&ATOMS' - CALL parser_init ( globenv % input_file_name, label, ierror, globenv ) - IF ( ierror /= 0 ) THEN - - IF( globenv % ionode) THEN - WRITE ( iw, '( A )' ) ' ATOM| No input section &ATOMS found ' - IF ( filen == ' ' ) THEN - CALL xstring ( globenv % project_name, ia, ie ) - filen = globenv % project_name ( ia:ie ) // '.dat' - END IF - ia = MIN ( LEN ( filen ), 20 ) - WRITE ( iw, '( A, T61, A )' ) ' ATOM| Try to read default file ', & - ADJUSTR ( filen ( 1:ia ) ) - CALL read_file ( filen, atype ) - END IF -!..broadcast the input data to all nodes -#if defined(__parallel) - source = globenv % source - allgrp = globenv % group - CALL mp_bcast ( atype % box, source, allgrp ) - CALL mp_bcast ( atype % ptype, source, allgrp ) - CALL mp_bcast ( atype % rtype, source, allgrp ) - CALL mp_bcast ( atype % n, source, allgrp ) - IF ( .NOT. globenv % ionode ) THEN - ALLOCATE ( atype % c(1:3,1:atype % n ), STAT = ios ) - IF ( ios /= 0 ) CALL stop_memory ( 'ATOM', 'atype%c', 3 * atype % n ) - END IF - CALL mp_bcast ( atype % c, source, allgrp ) - IF ( atype % rtype == 'POSVEL' ) THEN - IF ( .NOT. globenv % ionode ) THEN - ALLOCATE (atype % v(1:3,1:atype % n),STAT=ios) - IF ( ios /= 0 ) CALL stop_memory ( 'ATOM', 'atype%v', 3 * atype % n ) - END IF - CALL mp_bcast ( atype % v, source, allgrp ) - END IF -#endif - - ELSE - - CALL stop_prg ( 'atom','this part of the code not yet written' ) - - END IF - CALL parser_end -!..end of parsing the input section - -!..write some information to output - IF (globenv%ionode .and. globenv%print_level>0) THEN - WRITE ( iw, '( A )' ) ' ATOM| Box parameters [Angstrom]' - WRITE ( iw, '( A, T36, 3F15.5 )' ) & - ' ATOM| ', ( atype % box ( 1, i ), i = 1, 3 ) - WRITE ( iw, '( A, T36, 3F15.5 )' ) & - ' ATOM| ', ( atype % box ( 2, i ), i = 1, 3 ) - WRITE ( iw, '( A, T36, 3F15.5 )' ) & - ' ATOM| ', ( atype % box ( 3, i ), i = 1, 3 ) - WRITE ( iw, '( A, T71, I10 )' ) & - ' ATOM| Number of atoms read ', atype % n - IF ( globenv % print_level > 4 ) THEN - IF ( atype % rtype == 'POS' ) THEN - CALL print_c ( iw, atype % c ) - ELSE IF ( atype % rtype == 'POSVEL' ) THEN - CALL print_cv ( iw, atype % c, atype % v ) - END IF - END IF - WRITE ( iw, '()' ) - END IF - -END SUBROUTINE read_coord_vel - -!****************************************************************************** - -SUBROUTINE read_file ( filen, atype ) - - IMPLICIT NONE - -! Arguments - CHARACTER ( LEN = * ) :: filen - TYPE ( system_type ), INTENT ( INOUT ) :: atype - -! Locals - INTEGER :: iunit, i, j, ios - LOGICAL :: exists - -!------------------------------------------------------------------------------ - - INQUIRE ( FILE = filen, EXIST = exists ) - IF ( exists ) THEN - iunit = get_unit() - OPEN ( iunit, file = filen ) - READ ( iunit, * ) atype % n - - ALLOCATE ( atype % c ( 1:3, 1:atype % n ), STAT = ios ) - IF ( ios /= 0 ) CALL stop_memory ( 'ATOM', 'atype%c', 3 * atype % n ) - NULLIFY ( atype % v ) - - IF ( atype % rtype == 'POS' ) THEN - DO i = 1, atype % n - READ ( iunit, * ) atype % c ( 1:3, i ) - END DO - READ ( iunit, * ) atype % box ( 1, 1:3 ) - READ ( iunit, * ) atype % box ( 2, 1:3 ) - READ ( iunit, * ) atype % box ( 3, 1:3 ) - - ELSE IF ( atype % rtype == 'POSVEL' ) THEN - ALLOCATE ( atype % v ( 1:3, 1:atype % n ), STAT = ios ) - IF ( ios /= 0 ) & - CALL stop_memory ( 'ATOM', 'atype%v', 3 * atype % n ) - - DO i = 1, atype % n - READ ( iunit, * ) atype % c ( 1:3, i ) - END DO - READ ( iunit, * ) atype % box ( 1, 1:3 ) - READ ( iunit, * ) atype % box ( 2, 1:3 ) - READ ( iunit, * ) atype % box ( 3, 1:3 ) - DO i = 1, atype % n - READ ( iunit, * ) atype % v ( 1:3, i ) - END DO - - ELSE - CALL stop_prg ( 'ATOM', 'this rtype not programmed' ) - END IF - - ELSE - - CALL stop_prg ( 'ATOM', 'No information on atoms found ' ) - - END IF - -END SUBROUTINE read_file - -!****************************************************************************** - -SUBROUTINE print_c ( iw, c ) - - IMPLICIT NONE - -! Arguments - INTEGER, INTENT ( IN ) :: iw - REAL ( dbl ), DIMENSION ( :, : ), INTENT ( IN ) :: c - -! Locals - INTEGER :: i, n - -!------------------------------------------------------------------------------ - - n = SIZE ( c, 2 ) - WRITE ( iw, '( A )' ) ' ATOM| Atom coordinates [Angstrom]' - DO i = 1, n - WRITE ( iw, '( A, T26, I10, 3F15.5 )' ) ' ATOM| ', i, c ( :, i ) - END DO - -END SUBROUTINE print_c - -!****************************************************************************** - -SUBROUTINE print_cv ( iw, c, v ) - - IMPLICIT NONE - -! Arguments - INTEGER, INTENT ( IN ) :: iw - REAL ( dbl ), DIMENSION ( :, : ), INTENT ( IN ) :: c - REAL ( dbl ), DIMENSION ( :, : ), INTENT ( IN ) :: v - -! Locals - INTEGER :: i, n - -!------------------------------------------------------------------------------ - - n = SIZE ( c, 2 ) - WRITE ( iw, '( A )' ) ' ATOM| Atom coordinates [Angstrom]' - DO i = 1, n - WRITE ( iw, '( A, T8, I8, 3F10.4, 5X, 3F10.4 )' ) & - ' ATOM| ', i, c ( :, i ), v ( :, i ) - END DO - -END SUBROUTINE print_cv - -!****************************************************************************** - -END MODULE atoms_input diff --git a/src/atoms_input.cpp b/src/atoms_input.cpp new file mode 100644 index 0000000..e78e961 --- /dev/null +++ b/src/atoms_input.cpp @@ -0,0 +1,162 @@ +#include +#include + +typedef struct { + double box[3][3]; + char ptype[4]; + char rtype[21]; + int n; + double** c; + double** v; +} system_type; + +void read_coord_vel(system_type* atype, char* filen) { + int ierror, iw, source, allgrp, ia, ie, i, ios; + char label[6]; + + iw = globenv.scr; + + //..defaults + //..parse the input section + strcpy(label, "&ATOMS"); + parser_init(globenv.input_file_name, label, &ierror, globenv); + if (ierror != 0) { + if (globenv.ionode) { + printf("ATOM| No input section &ATOMS found\n"); + if (strcmp(filen, " ") == 0) { + ia = strlen(globenv.project_name); + ie = ia; + } else { + ia = strlen(filen); + ie = ia; + } + strcat(filen, ".dat"); + printf("ATOM| Try to read default file %s\n", filen); + read_file(filen, atype); + } + //..broadcast the input data to all nodes + #if defined(__parallel) + source = globenv.source; + allgrp = globenv.group; + mp_bcast(atype->box, source, allgrp); + mp_bcast(atype->ptype, source, allgrp); + mp_bcast(atype->rtype, source, allgrp); + mp_bcast(&atype->n, source, allgrp); + if (!globenv.ionode) { + atype->c = (double**)malloc(3 * sizeof(double*)); + for (i = 0; i < 3; i++) { + atype->c[i] = (double*)malloc(atype->n * sizeof(double)); + } + if (atype->c == NULL) { + stop_memory("ATOM", "atype%c", 3 * atype->n); + } + } + mp_bcast(atype->c, source, allgrp); + if (strcmp(atype->rtype, "POSVEL") == 0) { + if (!globenv.ionode) { + atype->v = (double**)malloc(3 * sizeof(double*)); + for (i = 0; i < 3; i++) { + atype->v[i] = (double*)malloc(atype->n * sizeof(double)); + } + if (atype->v == NULL) { + stop_memory("ATOM", "atype%v", 3 * atype->n); + } + } + mp_bcast(atype->v, source, allgrp); + } + #endif + } else { + stop_prg("atom", "this part of the code not yet written"); + } + parser_end(); + //..end of parsing the input section + + //..write some information to output + if (globenv.ionode && globenv.print_level > 0) { + printf("ATOM| Box parameters [Angstrom]\n"); + printf("ATOM| %15.5f %15.5f %15.5f\n", atype->box[0][0], atype->box[0][1], atype->box[0][2]); + printf("ATOM| %15.5f %15.5f %15.5f\n", atype->box[1][0], atype->box[1][1], atype->box[1][2]); + printf("ATOM| %15.5f %15.5f %15.5f\n", atype->box[2][0], atype->box[2][1], atype->box[2][2]); + printf("ATOM| Number of atoms read %d\n", atype->n); + if (globenv.print_level > 4) { + if (strcmp(atype->rtype, "POS") == 0) { + print_c(iw, atype->c); + } else if (strcmp(atype->rtype, "POSVEL") == 0) { + print_cv(iw, atype->c, atype->v); + } + } + printf("\n"); + } +} + +void read_file(char* filen, system_type* atype) { + int iunit, i, j, ios; + int exists; + + exists = access(filen, F_OK) != -1; + if (exists) { + iunit = get_unit(); + FILE* file = fopen(filen, "r"); + fscanf(file, "%d", &atype->n); + + atype->c = (double**)malloc(3 * sizeof(double*)); + for (i = 0; i < 3; i++) { + atype->c[i] = (double*)malloc(atype->n * sizeof(double)); + } + if (atype->c == NULL) { + stop_memory("ATOM", "atype%c", 3 * atype->n); + } + atype->v = NULL; + + if (strcmp(atype->rtype, "POS") == 0) { + for (i = 0; i < atype->n; i++) { + fscanf(file, "%lf %lf %lf", &atype->c[0][i], &atype->c[1][i], &atype->c[2][i]); + } + fscanf(file, "%lf %lf %lf", &atype->box[0][0], &atype->box[0][1], &atype->box[0][2]); + fscanf(file, "%lf %lf %lf", &atype->box[1][0], &atype->box[1][1], &atype->box[1][2]); + fscanf(file, "%lf %lf %lf", &atype->box[2][0], &atype->box[2][1], &atype->box[2][2]); + } else if (strcmp(atype->rtype, "POSVEL") == 0) { + atype->v = (double**)malloc(3 * sizeof(double*)); + for (i = 0; i < 3; i++) { + atype->v[i] = (double*)malloc(atype->n * sizeof(double)); + } + if (atype->v == NULL) { + stop_memory("ATOM", "atype%v", 3 * atype->n); + } + + for (i = 0; i < atype->n; i++) { + fscanf(file, "%lf %lf %lf", &atype->c[0][i], &atype->c[1][i], &atype->c[2][i]); + } + fscanf(file, "%lf %lf %lf", &atype->box[0][0], &atype->box[0][1], &atype->box[0][2]); + fscanf(file, "%lf %lf %lf", &atype->box[1][0], &atype->box[1][1], &atype->box[1][2]); + fscanf(file, "%lf %lf %lf", &atype->box[2][0], &atype->box[2][1], &atype->box[2][2]); + for (i = 0; i < atype->n; i++) { + fscanf(file, "%lf %lf %lf", &atype->v[0][i], &atype->v[1][i], &atype->v[2][i]); + } + } else { + stop_prg("ATOM", "this rtype not programmed"); + } + + fclose(file); + } else { + stop_prg("ATOM", "No information on atoms found"); + } +} + +void print_c(int iw, double** c, int n) { + int i; + + printf("ATOM| Atom coordinates [Angstrom]\n"); + for (i = 0; i < n; i++) { + printf("ATOM| %d %15.5f %15.5f %15.5f\n", i, c[0][i], c[1][i], c[2][i]); + } +} + +void print_cv(int iw, double** c, double** v, int n) { + int i; + + printf("ATOM| Atom coordinates [Angstrom]\n"); + for (i = 0; i < n; i++) { + printf("ATOM| %d %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n", i, c[0][i], c[1][i], c[2][i], v[0][i], v[1][i], v[2][i]); + } +} diff --git a/src/header.F b/src/header.F deleted file mode 100644 index af588a6..0000000 --- a/src/header.F +++ /dev/null @@ -1,159 +0,0 @@ -!-----------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright (C) 2000 CP2K developers group ! -!-----------------------------------------------------------------------------! - -MODULE header - - IMPLICIT NONE - - PRIVATE - PUBLIC :: fist_header, tbmd_header, qs_header, wave_header, faust_header - -CONTAINS - -!****************************************************************************** - -SUBROUTINE fist_header ( iw ) - - IMPLICIT NONE - -! Arguments - INTEGER :: iw - -!------------------------------------------------------------------------------ - - WRITE ( iw, '( / )' ) - WRITE ( iw, '( 11(10x,a,/) )' ) & - ' _ __ ', & - ' ___/ \\/ \\__ ', & - ' /| | | | \\ ', & - ' / | | | | | ', & - ' | | | | | | ', & - ' FRONTIERS IN | \\__| |__|__| SIMULATION TECHNOLOGY', & - ' | \\__/ \\ ', & - ' C.J. Mundy \\ \\____/ | ', & - ' S. Balasubramanian / | 1998-1999', & - ' Ken Bagchi \\ / Version 0.0', & - ' ' - -END SUBROUTINE fist_header - -!****************************************************************************** - -SUBROUTINE tbmd_header ( iw ) - IMPLICIT NONE - INTEGER :: iw - -!------------------------------------------------------------------------------ - - WRITE ( iw, '( / )' ) - WRITE ( iw, '( 14(12x,a,/) )' ) & - ' ********************** ************* ******** ', & - ' ************************ *************** ********** ', & - ' *************** *** *************** *********** ', & - ' **** **** **** **** **** *** **** ', & - ' **** *********** **** **** *** **** ', & - ' **** **** *** **** **** *********** ', & - ' **** *********** **** **** ********** ', & - ' **** ********* **** **** ******* ', & - ' ', & - ' University of Zurich ', & - ' 2000 ', & - ' ', & - ' Version 0.0 ', & - ' ' - -END SUBROUTINE tbmd_header - -!****************************************************************************** - -SUBROUTINE qs_header ( iw ) - - IMPLICIT NONE - -! Arguments - INTEGER :: iw - - WRITE ( iw, '( / )' ) - WRITE ( iw, '( 14(5x,a,/) )' ) & - ' ***** ',& - ' ********* ** *** ** ',& - ' **** **** *** *** ****** ',& - '**** **** ** *** ***** *** *** *********** ****** ****** ',& - '**** ******** *** *** ******* ****** **** **** *** *** ** ***',& - ' **** ***** *** *** *** *** **** ****** **** ******** *******',& - ' ********** ******* *** ******* ****** ******** **** ***** ',& - ' ***** ** ***** *** ***** *** ******* **** ****** *** ',& - ' *** ',& - ' ',& - ' MPI Festkoerperforschung Stuttgart ',& - ' 1999 ',& - ' ',& - ' Version 0.0 ',& - ' ' -END SUBROUTINE qs_header - -!****************************************************************************** - -SUBROUTINE wave_header ( iw ) - - IMPLICIT NONE - -! Arguments - INTEGER :: iw - -!------------------------------------------------------------------------------ - - WRITE ( iw, '( / )' ) - WRITE ( iw, '( 14(12x,a,/) )' ) & - ' **** **** ********* **** *********** ', & - ' **** **** *********** **** ************* ', & - ' **** *** ******* ******** ******** ', & - ' **** ***** ******* ******** **** ************ ', & - ' **** *** *** ******************* **** ************ ', & - ' ******** **************************** **** ', & - ' ****** ****** **** **** ******* ********** ', & - ' **** **** **** **** ****** ******** ', & - ' ', & - ' MPI Festkoerperforschung Stuttgart ', & - ' 1999 ', & - ' ', & - ' Version 0.0 ', & - ' ' - -END SUBROUTINE wave_header - -!****************************************************************************** - -SUBROUTINE faust_header ( iw ) - - IMPLICIT NONE - -! Arguments - INTEGER :: iw - -!------------------------------------------------------------------------------ - - WRITE (iw,'(/)') - WRITE (iw,'(14(12x,a,/))') & - ' **** **** ********* **** *********** ', & - ' **** **** *********** **** ************* ', & - ' **** *** ******* ******** ******** ', & - ' **** ***** ******* ******** **** ************ ', & - ' **** *** *** ******************* **** ************ ', & - ' ******** **************************** **** ', & - ' ****** ****** **** **** ******* ********** ', & - ' **** **** **** **** ****** ******** ', & - ' ', & - ' MPI Festkoerperforschung Stuttgart ', & - ' 1999 ', & - ' ', & - ' Version 0.0 ', & - ' ' - -END SUBROUTINE faust_header - -!****************************************************************************** - -END MODULE header diff --git a/src/header.c b/src/header.c new file mode 100644 index 0000000..719853f --- /dev/null +++ b/src/header.c @@ -0,0 +1,119 @@ +/*---------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright (C) 2000 CP2K developers group */ +/*---------------------------------------------------------------------------*/ + +#include "header.h" + +//***************************************************************************** + +void fist_header(FILE *iw) { + +//----------------------------------------------------------------------------- + + fprintf(iw, "\n"); + fprintf(iw, " _ __ \n"); + fprintf(iw, " ___/ \\/ \\__ \n"); + fprintf(iw, " /| | | | \\ \n"); + fprintf(iw, " / | | | | | \n"); + fprintf(iw, " | | | | | | \n"); + fprintf(iw, " FRONTIERS IN | \\__| |__|__| SIMULATION TECHNOLOGY\n"); + fprintf(iw, " | \\__/ \\ \n"); + fprintf(iw, " C.J. Mundy \\ \\____/ | \n"); + fprintf(iw, " S. Balasubramanian / | 1998-1999\n"); + fprintf(iw, " Ken Bagchi \\ / Version 0.0\n"); + fprintf(iw, " \n"); +} + +//***************************************************************************** + +void tbmd_header(FILE *iw) { + +//----------------------------------------------------------------------------- + + fprintf(iw, "\n"); + fprintf(iw, " ********************** ************* ******** \n"); + fprintf(iw, " ************************ *************** ********** \n"); + fprintf(iw, " *************** *** *************** *********** \n"); + fprintf(iw, " **** **** **** **** **** *** **** \n"); + fprintf(iw, " **** *********** **** **** *** **** \n"); + fprintf(iw, " **** **** *** **** **** *********** \n"); + fprintf(iw, " **** *********** **** **** ********** \n"); + fprintf(iw, " **** ********* **** **** ******* \n"); + fprintf(iw, " \n"); + fprintf(iw, " University of Zurich \n"); + fprintf(iw, " 2000 \n"); + fprintf(iw, " \n"); + fprintf(iw, " Version 0.0 \n"); + fprintf(iw, " \n"); +} + +//***************************************************************************** + +void qs_header(FILE *iw) { + + fprintf(iw, "\n"); + fprintf(iw, " ***** \n"); + fprintf(iw, " ********* ** *** ** \n"); + fprintf(iw, " **** **** *** *** ****** \n"); + fprintf(iw, "**** **** ** *** ***** *** *** *********** ****** ****** \n"); + fprintf(iw, "**** ******** *** *** ******* ****** **** **** *** *** ** ***\n"); + fprintf(iw, " **** ***** *** *** *** *** **** ****** **** ******** *******\n"); + fprintf(iw, " ********** ******* *** ******* ****** ******** **** ***** \n"); + fprintf(iw, " ***** ** ***** *** ***** *** ******* **** ****** *** \n"); + fprintf(iw, " *** \n"); + fprintf(iw, " \n"); + fprintf(iw, " MPI Festkoerperforschung Stuttgart \n"); + fprintf(iw, " 1999 \n"); + fprintf(iw, " \n"); + fprintf(iw, " Version 0.0 \n"); + fprintf(iw, " \n"); +} + +//***************************************************************************** + +void wave_header(FILE *iw) { + +//----------------------------------------------------------------------------- + + fprintf(iw, "\n"); + fprintf(iw, " **** **** ********* **** *********** \n"); + fprintf(iw, " **** **** *********** **** ************* \n"); + fprintf(iw, " **** *** ******* ******** ******** \n"); + fprintf(iw, " **** ***** ******* ******** **** ************ \n"); + fprintf(iw, " **** *** *** ******************* **** ************ \n"); + fprintf(iw, " ******** **************************** **** \n"); + fprintf(iw, " ****** ****** **** **** ******* ********** \n"); + fprintf(iw, " **** **** **** **** ****** ******** \n"); + fprintf(iw, " \n"); + fprintf(iw, " MPI Festkoerperforschung Stuttgart \n"); + fprintf(iw, " 1999 \n"); + fprintf(iw, " \n"); + fprintf(iw, " Version 0.0 \n"); + fprintf(iw, " \n"); +} + +//***************************************************************************** + +void faust_header(FILE *iw) { + +//----------------------------------------------------------------------------- + + fprintf(iw, "\n"); + fprintf(iw, " **** **** ********* **** *********** \n"); + fprintf(iw, " **** **** *********** **** ************* \n"); + fprintf(iw, " **** *** ******* ******** ******** \n"); + fprintf(iw, " **** ***** ******* ******** **** ************ \n"); + fprintf(iw, " **** *** *** ******************* **** ************ \n"); + fprintf(iw, " ******** **************************** **** \n"); + fprintf(iw, " ****** ****** **** **** ******* ********** \n"); + fprintf(iw, " **** **** **** **** ****** ******** \n"); + fprintf(iw, " \n"); + fprintf(iw, " MPI Festkoerperforschung Stuttgart \n"); + fprintf(iw, " 1999 \n"); + fprintf(iw, " \n"); + fprintf(iw, " Version 0.0 \n"); + fprintf(iw, " \n"); +} + +/*****************************************************************************/ diff --git a/src/header.h b/src/header.h new file mode 100644 index 0000000..82c10f4 --- /dev/null +++ b/src/header.h @@ -0,0 +1,21 @@ +#ifndef _HEADER_H +#define _HEADER_H + +// Secure block +#ifdef __cplusplus +extern "C" { +#endif + +#include + +void fist_header(FILE *iw); +void tbmd_header(FILE *iw); +void qs_header(FILE *iw); +void wave_header(FILE *iw); +void faust_header(FILE *iw); + +#ifdef __cplusplus +} +#endif + +#endif \ No newline at end of file diff --git a/src/kinds.cpp b/src/kinds.cpp index 3dcd63a..5106c4c 100644 --- a/src/kinds.cpp +++ b/src/kinds.cpp @@ -13,6 +13,7 @@ - HITACHI */ +#include #include namespace kinds { diff --git a/src/slater_koster.h b/src/slater_koster.h new file mode 100644 index 0000000..4217ee4 --- /dev/null +++ b/src/slater_koster.h @@ -0,0 +1,36 @@ +#ifndef _SLATER_KOSTER_H +#define _SLATER_KOSTER_H + +#ifdef __cplusplus +extern "C" { +#endif + +const double zero = 0.0; +const double one = 1.0; + +void sph(int l, double r[3], double gsl[7]); + +void dsph(int l, double r[3], double dgsl[7][3]); + +void out_prod(double mat[][3], double v1[], double v2[]); + +void out_dprod(double dmat[][3], double v1[], double v2[], double dv1[], + double dv2[]); + +double dpro(double g1[][7], double g2[][7], int l1, int l2); + +void m_integrals(); + +void gmat(int l1, int l2, double r[3], double gmu[7][7][4], double dll[7][7][7]); + +void dgmat(int l1, int l2, double r[3], double dgmu[7][7][4][3], double ddll[7][7][7][3]); + +void sk_test(); + +void dsk_test(); + +#ifdef __cplusplus +} +#endif + +#endif \ No newline at end of file diff --git a/src/slater_koster_matr.F b/src/slater_koster_matr.F deleted file mode 100644 index f780153..0000000 --- a/src/slater_koster_matr.F +++ /dev/null @@ -1,1509 +0,0 @@ -!------------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright (C) 2000 CP2K developers group ! -!------------------------------------------------------------------------------! -!!> Calculation of two-center s-f Slater-Koster integrals -!! A.K. McMahan, Phys. Rev. B, 58, p4293 (1998) -!! -!! Sign convention: t(lp,lq)=(-1)**(lp+lq)*(lp m,lq m) -!! where (lp m,lq m) is the integral in the diatomic coordinate system -!! -!!< - MODULE slater_koster_matr -!------------------------------------------------------------------------------! - USE kinds, ONLY : dbl - USE slater_koster_util, ONLY : sph, dsph, out_prod, out_dprod -! - IMPLICIT NONE -! - PRIVATE - PUBLIC :: gmat, dgmat -! - REAL (dbl), PARAMETER :: zero = 0._dbl - REAL (dbl), PARAMETER :: one = 1._dbl -!------------------------------------------------------------------------------! -! - CONTAINS -! -!------------------------------------------------------------------------------! - SUBROUTINE gmat(l1,l2,r,gmu,dll) - IMPLICIT NONE - INTEGER, INTENT (IN) :: l1, l2 - REAL (dbl), INTENT (IN) :: r(3) - REAL (dbl), INTENT (OUT) :: gmu(0:6,0:6,0:3) - REAL (dbl), INTENT (OUT) :: dll(0:6,0:6,0:6) - - REAL (dbl) :: gsl(0:6,0:3) - INTEGER :: mu, mo, m1, m2, i - - mu = min(l1,l2) - mo = max(l1,l2) - m1 = 2*l1 - m2 = 2*l2 - gsl(0:6,0:3) = zero - DO i = 0, mo - CALL sph(i,r,gsl(0:6,i)) - END DO - dll(0:6,0:6,0) = zero - DO i = 0, 2*mo - dll(i,i,0) = 1._dbl - END DO - - IF (l1==0) THEN -!s-l - gmu(0,0:m2,0) = gsl(0:m2,l2) - ELSE IF (l2==0) THEN -!l-s - gmu(0:m1,0,0) = gsl(0:m1,l1) - ELSE - SELECT CASE (l1) - CASE (1) - SELECT CASE (l2) - CASE (1) -!p-p - CALL out_prod(dll(0:2,0:2,2),gsl(0:2,1),gsl(0:2,1)) - gmu(0:2,0:2,0) = dll(0:2,0:2,2) - gmu(0:2,0:2,1) = -dll(0:2,0:2,2) + dll(0:2,0:2,0) - CASE (2) -!p-d - CALL out_prod(dll(0:2,0:4,3),gsl(0:2,1),gsl(0:4,2)) - CALL set_mat(dll(0:2,0:4,1),r,1) - gmu(0:2,0:4,0) = dll(0:2,0:4,3) - gmu(0:2,0:4,1) = -2._dbl/sqrt(3._dbl)*dll(0:2,0:4,3) + & - dll(0:2,0:4,1) - CASE (3) -!p-f - CALL out_prod(dll(0:2,0:6,4),gsl(0:2,1),gsl(0:6,3)) - CALL set_mat(dll(0:2,0:6,2),r,2) - gmu(0:2,0:6,0) = dll(0:2,0:6,4) - gmu(0:2,0:6,1) = -sqrt(1.5_dbl)*dll(0:2,0:6,4) + & - 0.5_dbl*dll(0:2,0:6,2) - END SELECT - CASE (2) - SELECT CASE (l2) - CASE (1) -!d-p - CALL out_prod(dll(0:4,0:2,3),gsl(0:4,2),gsl(0:2,1)) - CALL set_mat(dll(0:4,0:2,1),r,1) - gmu(0:4,0:2,0) = dll(0:4,0:2,3) - gmu(0:4,0:2,1) = -2._dbl/sqrt(3._dbl)*dll(0:4,0:2,3) + & - dll(0:4,0:2,1) - CASE (2) -!d-d - CALL out_prod(dll(0:4,0:4,4),gsl(0:4,2),gsl(0:4,2)) - CALL set_mat(dll(0:4,0:4,2),r,2) - gmu(0:4,0:4,0) = dll(0:4,0:4,4) - gmu(0:4,0:4,1) = -4._dbl/3._dbl*dll(0:4,0:4,4) + & - dll(0:4,0:4,2) + dll(0:4,0:4,0) - gmu(0:4,0:4,2) = 1._dbl/3._dbl*dll(0:4,0:4,4) - dll(0:4,0:4,2) - CASE (3) -!d-f - CALL out_prod(dll(0:4,0:6,5),gsl(0:4,2),gsl(0:6,3)) - CALL set_mat(dll(0:4,0:6,3),r,3) - CALL set_mat(dll(0:4,0:6,1),r,1) - gmu(0:4,0:6,0) = dll(0:4,0:6,5) - gmu(0:4,0:6,1) = -sqrt(2._dbl)*dll(0:4,0:6,5) + & - 0.5_dbl*dll(0:4,0:6,3) + 0.5_dbl*dll(0:4,0:6,1) - gmu(0:4,0:6,2) = sqrt(0.2_dbl)*dll(0:4,0:6,5) - & - sqrt(0.4_dbl)*dll(0:4,0:6,3) - END SELECT - CASE (3) - SELECT CASE (l2) - CASE (1) -!f-p - CALL out_prod(dll(0:6,0:2,4),gsl(0:6,3),gsl(0:2,1)) - CALL set_mat(dll(0:6,0:2,2),r,2) - gmu(0:6,0:2,0) = dll(0:6,0:2,4) - gmu(0:6,0:2,1) = -sqrt(1.5_dbl)*dll(0:6,0:2,4) + & - 0.5_dbl*dll(0:6,0:2,2) - CASE (2) -!f-d - CALL out_prod(dll(0:6,0:4,5),gsl(0:6,3),gsl(0:4,2)) - CALL set_mat(dll(0:6,0:4,3),r,3) - CALL set_mat(dll(0:6,0:4,1),r,1) - gmu(0:6,0:4,0) = dll(0:6,0:4,5) - gmu(0:6,0:4,1) = -sqrt(2._dbl)*dll(0:6,0:4,5) + & - 0.5_dbl*dll(0:6,0:4,3) + 0.5_dbl*dll(0:6,0:4,1) - gmu(0:6,0:4,2) = sqrt(0.2_dbl)*dll(0:6,0:4,5) - & - sqrt(0.4_dbl)*dll(0:6,0:4,3) - CASE (3) -!f-f - CALL out_prod(dll(0:6,0:6,6),gsl(0:6,3),gsl(0:6,3)) - CALL set_mat(dll(0:6,0:6,4),r,4) - CALL set_mat(dll(0:6,0:6,2),r,2) - gmu(0:6,0:6,0) = dll(0:6,0:6,6) - gmu(0:6,0:6,1) = -1.5_dbl*dll(0:6,0:6,6) + & - 0.625_dbl*dll(0:6,0:6,4) + 0.625_dbl*dll(0:6,0:6,2) + & - 0.625_dbl*dll(0:6,0:6,0) - gmu(0:6,0:6,2) = 0.6_dbl*dll(0:6,0:6,6) - dll(0:6,0:6,4) - gmu(0:6,0:6,3) = -0.1_dbl*dll(0:6,0:6,6) + & - 0.375_dbl*dll(0:6,0:6,4) - 0.625_dbl*dll(0:6,0:6,2) + & - 0.375_dbl*dll(0:6,0:6,0) - END SELECT - END SELECT - END IF - - END SUBROUTINE gmat -!------------------------------------------------------------------------------! - SUBROUTINE dgmat(l1,l2,r,dgmu,ddll) - IMPLICIT NONE - INTEGER, INTENT (IN) :: l1, l2 - REAL (dbl), INTENT (IN) :: r(3) - REAL (dbl), INTENT (OUT) :: dgmu(0:6,0:6,0:3,1:3) - REAL (dbl), INTENT (OUT) :: ddll(0:6,0:6,0:6,1:3) - - REAL (dbl) :: gsl(0:6,0:3) - REAL (dbl) :: dgsl(0:6,0:3,1:3) - INTEGER :: mu, mo, m1, m2, i - - mu = min(l1,l2) - mo = max(l1,l2) - m1 = 2*l1 - m2 = 2*l2 - gsl(0:6,0:3) = zero - dgsl(0:6,0:3,1:3) = zero - DO i = 0, mo - CALL sph(i,r,gsl(0:6,i)) - CALL dsph(i,r,dgsl(0:6,i,1:3)) - END DO - - IF (l1==0) THEN -!s-l - dgmu(0,0:m2,0,1:3) = dgsl(0:m2,l2,1:3) - ELSE IF (l2==0) THEN -!l-s - dgmu(0:m1,0,0,1:3) = dgsl(0:m1,l1,1:3) - ELSE - SELECT CASE (l1) - CASE (1) - SELECT CASE (l2) - CASE (1) -!p-p - CALL out_dprod(ddll(0:2,0:2,2,1),gsl(0:2,1),gsl(0:2,1), & - dgsl(0:2,1,1),dgsl(0:2,1,1)) - CALL out_dprod(ddll(0:2,0:2,2,2),gsl(0:2,1),gsl(0:2,1), & - dgsl(0:2,1,2),dgsl(0:2,1,2)) - CALL out_dprod(ddll(0:2,0:2,2,3),gsl(0:2,1),gsl(0:2,1), & - dgsl(0:2,1,3),dgsl(0:2,1,3)) - dgmu(0:2,0:2,0,1:3) = ddll(0:2,0:2,2,1:3) - dgmu(0:2,0:2,1,1:3) = -ddll(0:2,0:2,2,1:3) - CASE (2) -!p-d - CALL out_dprod(ddll(0:2,0:4,3,1),gsl(0:2,1),gsl(0:4,2), & - dgsl(0:2,1,1),dgsl(0:4,2,1)) - CALL out_dprod(ddll(0:2,0:4,3,2),gsl(0:2,1),gsl(0:4,2), & - dgsl(0:2,1,2),dgsl(0:4,2,2)) - CALL out_dprod(ddll(0:2,0:4,3,3),gsl(0:2,1),gsl(0:4,2), & - dgsl(0:2,1,3),dgsl(0:4,2,3)) - CALL set_dmat(ddll(0:2,0:4,1,1:3),r,1) - dgmu(0:2,0:4,0,1:3) = ddll(0:2,0:4,3,1:3) - dgmu(0:2,0:4,1,1:3) = -2._dbl/sqrt(3._dbl)*ddll(0:2,0:4,3,1:3) + & - ddll(0:2,0:4,1,1:3) - CASE (3) -!p-f - CALL out_dprod(ddll(0:2,0:6,4,1),gsl(0:2,1),gsl(0:6,3), & - dgsl(0:2,1,1),dgsl(0:6,3,1)) - CALL out_dprod(ddll(0:2,0:6,4,2),gsl(0:2,1),gsl(0:6,3), & - dgsl(0:2,1,2),dgsl(0:6,3,2)) - CALL out_dprod(ddll(0:2,0:6,4,3),gsl(0:2,1),gsl(0:6,3), & - dgsl(0:2,1,3),dgsl(0:6,3,3)) - CALL set_dmat(ddll(0:2,0:6,2,1:3),r,2) - dgmu(0:2,0:6,0,1:3) = ddll(0:2,0:6,4,1:3) - dgmu(0:2,0:6,1,1:3) = -sqrt(1.5_dbl)*ddll(0:2,0:6,4,1:3) + & - 0.5_dbl*ddll(0:2,0:6,2,1:3) - END SELECT - CASE (2) - SELECT CASE (l2) - CASE (1) -!d-p - CALL out_dprod(ddll(0:4,0:2,3,1),gsl(0:4,2),gsl(0:2,1), & - dgsl(0:4,2,1),dgsl(0:2,1,1)) - CALL out_dprod(ddll(0:4,0:2,3,2),gsl(0:4,2),gsl(0:2,1), & - dgsl(0:4,2,2),dgsl(0:2,1,2)) - CALL out_dprod(ddll(0:4,0:2,3,3),gsl(0:4,2),gsl(0:2,1), & - dgsl(0:4,2,3),dgsl(0:2,1,3)) - CALL set_dmat(ddll(0:4,0:2,1,1:3),r,1) - dgmu(0:4,0:2,0,1:3) = ddll(0:4,0:2,3,1:3) - dgmu(0:4,0:2,1,1:3) = -2._dbl/sqrt(3._dbl)*ddll(0:4,0:2,3,1:3) + & - ddll(0:4,0:2,1,1:3) - CASE (2) -!d-d - CALL out_dprod(ddll(0:4,0:4,4,1),gsl(0:4,2),gsl(0:4,2), & - dgsl(0:4,2,1),dgsl(0:4,2,1)) - CALL out_dprod(ddll(0:4,0:4,4,2),gsl(0:4,2),gsl(0:4,2), & - dgsl(0:4,2,2),dgsl(0:4,2,2)) - CALL out_dprod(ddll(0:4,0:4,4,3),gsl(0:4,2),gsl(0:4,2), & - dgsl(0:4,2,3),dgsl(0:4,2,3)) - CALL set_dmat(ddll(0:4,0:4,2,1:3),r,2) - dgmu(0:4,0:4,0,1:3) = ddll(0:4,0:4,4,1:3) - dgmu(0:4,0:4,1,1:3) = -4._dbl/3._dbl*ddll(0:4,0:4,4,1:3) + & - ddll(0:4,0:4,2,1:3) - dgmu(0:4,0:4,2,1:3) = 1._dbl/3._dbl*ddll(0:4,0:4,4,1:3) - & - ddll(0:4,0:4,2,1:3) - CASE (3) -!d-f - CALL out_dprod(ddll(0:4,0:6,5,1),gsl(0:4,2),gsl(0:6,3), & - dgsl(0:4,2,1),dgsl(0:6,3,1)) - CALL out_dprod(ddll(0:4,0:6,5,2),gsl(0:4,2),gsl(0:6,3), & - dgsl(0:4,2,2),dgsl(0:6,3,2)) - CALL out_dprod(ddll(0:4,0:6,5,3),gsl(0:4,2),gsl(0:6,3), & - dgsl(0:4,2,3),dgsl(0:6,3,3)) - CALL set_dmat(ddll(0:4,0:6,3,1:3),r,3) - CALL set_dmat(ddll(0:4,0:6,1,1:3),r,1) - dgmu(0:4,0:6,0,1:3) = ddll(0:4,0:6,5,1:3) - dgmu(0:4,0:6,1,1:3) = -sqrt(2._dbl)*ddll(0:4,0:6,5,1:3) + & - 0.5_dbl*ddll(0:4,0:6,3,1:3) + 0.5_dbl*ddll(0:4,0:6,1,1:3) - dgmu(0:4,0:6,2,1:3) = sqrt(0.2_dbl)*ddll(0:4,0:6,5,1:3) - & - sqrt(0.4_dbl)*ddll(0:4,0:6,3,1:3) - END SELECT - CASE (3) - SELECT CASE (l2) - CASE (1) -!f-p - CALL out_dprod(ddll(0:6,0:2,4,1),gsl(0:6,3),gsl(0:2,1), & - dgsl(0:6,3,1),dgsl(0:2,1,1)) - CALL out_dprod(ddll(0:6,0:2,4,2),gsl(0:6,3),gsl(0:2,1), & - dgsl(0:6,3,2),dgsl(0:2,1,2)) - CALL out_dprod(ddll(0:6,0:2,4,3),gsl(0:6,3),gsl(0:2,1), & - dgsl(0:6,3,3),dgsl(0:2,1,3)) - CALL set_dmat(ddll(0:6,0:2,2,1:3),r,2) - dgmu(0:6,0:2,0,1:3) = ddll(0:6,0:2,4,1:3) - dgmu(0:6,0:2,1,1:3) = -sqrt(1.5_dbl)*ddll(0:6,0:2,4,1:3) + & - 0.5_dbl*ddll(0:6,0:2,2,1:3) - CASE (2) -!f-d - CALL out_dprod(ddll(0:6,0:4,5,1),gsl(0:6,3),gsl(0:4,2), & - dgsl(0:6,3,1),dgsl(0:4,2,1)) - CALL out_dprod(ddll(0:6,0:4,5,2),gsl(0:6,3),gsl(0:4,2), & - dgsl(0:6,3,2),dgsl(0:4,2,2)) - CALL out_dprod(ddll(0:6,0:4,5,3),gsl(0:6,3),gsl(0:4,2), & - dgsl(0:6,3,3),dgsl(0:4,2,3)) - CALL set_dmat(ddll(0:6,0:4,3,1:3),r,3) - CALL set_dmat(ddll(0:6,0:4,1,1:3),r,1) - dgmu(0:6,0:4,0,1:3) = ddll(0:6,0:4,5,1:3) - dgmu(0:6,0:4,1,1:3) = -sqrt(2._dbl)*ddll(0:6,0:4,5,1:3) + & - 0.5_dbl*ddll(0:6,0:4,3,1:3) + 0.5_dbl*ddll(0:6,0:4,1,1:3) - dgmu(0:6,0:4,2,1:3) = sqrt(0.2_dbl)*ddll(0:6,0:4,5,1:3) - & - sqrt(0.4_dbl)*ddll(0:6,0:4,3,1:3) - CASE (3) -!f-f - CALL out_dprod(ddll(0:6,0:6,6,1),gsl(0:6,3),gsl(0:6,3), & - dgsl(0:6,3,1),dgsl(0:6,3,1)) - CALL out_dprod(ddll(0:6,0:6,6,2),gsl(0:6,3),gsl(0:6,3), & - dgsl(0:6,3,2),dgsl(0:6,3,2)) - CALL out_dprod(ddll(0:6,0:6,6,3),gsl(0:6,3),gsl(0:6,3), & - dgsl(0:6,3,3),dgsl(0:6,3,3)) - CALL set_dmat(ddll(0:6,0:6,4,1:3),r,4) - CALL set_dmat(ddll(0:6,0:6,2,1:3),r,2) - dgmu(0:6,0:6,0,1:3) = ddll(0:6,0:6,6,1:3) - dgmu(0:6,0:6,1,1:3) = -1.5_dbl*ddll(0:6,0:6,6,1:3) + & - 0.625_dbl*ddll(0:6,0:6,4,1:3) + 0.625_dbl*ddll(0:6,0:6,2,1:3) - dgmu(0:6,0:6,2,1:3) = 0.6_dbl*ddll(0:6,0:6,6,1:3) - ddll(0:6,0:6,4,1:3) - dgmu(0:6,0:6,3,1:3) = -0.1_dbl*ddll(0:6,0:6,6,1:3) + & - 0.375_dbl*ddll(0:6,0:6,4,1:3) - 0.625_dbl*ddll(0:6,0:6,2,1:3) - END SELECT - END SELECT - END IF - - END SUBROUTINE dgmat -!------------------------------------------------------------------------------! - SUBROUTINE set_mat(mat,r,num) - IMPLICIT NONE - REAL (dbl), INTENT (IN) :: r(3) - INTEGER, INTENT (IN) :: num - REAL (dbl), INTENT (OUT) :: mat(0:,0:) - INTEGER :: l1, l2 - REAL (dbl), PARAMETER :: s3 = 0.577350269189626_dbl ! sqrt(1._dbl/3._dbl) - REAL (dbl) :: x, y, z, xx, yy, zz, xy, xz, yz - - l1 = size(mat(:,0)) - l2 = size(mat(0,:)) - x = r(1) - y = r(2) - z = r(3) - IF (l1==5 .AND. l2==3 .AND. num==1) THEN -! d1(d,p) - mat(0:4,0:2) = zero - mat(0,0) = 2._dbl*s3*z - mat(1,0) = x - mat(2,0) = y - mat(0,1) = -s3*x - mat(1,1) = z - mat(3,1) = x - mat(4,1) = y - mat(0,2) = -s3*y - mat(2,2) = z - mat(3,2) = -y - mat(4,2) = x - ELSE IF (l1==3 .AND. l2==5 .AND. num==1) THEN -! d1(p,d) - mat(0:2,0:4) = zero - mat(0,0) = 2._dbl*s3*z - mat(0,1) = x - mat(0,2) = y - mat(1,0) = -s3*x - mat(1,1) = z - mat(1,3) = x - mat(1,4) = y - mat(2,0) = -s3*y - mat(2,2) = z - mat(2,3) = -y - mat(2,4) = x - ELSE - xx = x*x - yy = y*y - zz = z*z - xy = x*y - xz = x*z - yz = y*z - IF (l1==5 .AND. l2==5 .AND. num==2) THEN -! d2(d,d) - mat(0,0) = zz - 2._dbl/3._dbl - mat(0,1) = s3*xz - mat(0,2) = s3*yz - mat(0,3) = -s3*(xx-yy) - mat(0,4) = -2._dbl*s3*xy - mat(1,0) = mat(0,1) - mat(2,0) = mat(0,2) - mat(3,0) = mat(0,3) - mat(4,0) = mat(0,4) - mat(1,1) = -yy - mat(1,2) = xy - mat(1,3) = xz - mat(1,4) = yz - mat(2,1) = mat(1,2) - mat(3,1) = mat(1,3) - mat(4,1) = mat(1,4) - mat(2,2) = -xx - mat(2,3) = -yz - mat(2,4) = xz - mat(3,2) = mat(2,3) - mat(4,2) = mat(2,4) - mat(3,3) = -zz - mat(3,4) = zero - mat(4,3) = zero - mat(4,4) = -zz - ELSE IF (l1==7 .AND. l2==3 .AND. num==2) THEN -! d2(f,p) - mat(0,0) = sqrt(1.5_dbl)*(3._dbl*zz-1._dbl) - mat(0,1) = -sqrt(6._dbl)*xz - mat(0,2) = -sqrt(6._dbl)*yz - mat(1,0) = 4._dbl*xz - mat(1,1) = 0.5_dbl*(5._dbl*zz-2._dbl*xx-1._dbl) - mat(1,2) = -xy - mat(2,0) = 4._dbl*yz - mat(2,1) = -xy - mat(2,2) = 0.5_dbl*(5._dbl*zz-2._dbl*yy-1._dbl) - mat(3,0) = sqrt(2.5_dbl)*(xx-yy) - mat(3,1) = sqrt(10._dbl)*xz - mat(3,2) = -sqrt(10._dbl)*yz - mat(4,0) = sqrt(10._dbl)*xy - mat(4,1) = sqrt(10._dbl)*yz - mat(4,2) = sqrt(10._dbl)*xz - mat(5,0) = zero - mat(5,1) = sqrt(3.75_dbl)*(xx-yy) - mat(5,2) = -sqrt(15._dbl)*xy - mat(6,0) = zero - mat(6,1) = sqrt(15._dbl)*xy - mat(6,2) = sqrt(3.75_dbl)*(xx-yy) - ELSE IF (l1==3 .AND. l2==7 .AND. num==2) THEN -! d2(p,f) - mat(0,0) = sqrt(1.5_dbl)*(3._dbl*zz-1._dbl) - mat(1,0) = -sqrt(6._dbl)*xz - mat(2,0) = -sqrt(6._dbl)*yz - mat(0,1) = 4._dbl*xz - mat(1,1) = 0.5_dbl*(5._dbl*zz-2._dbl*xx-1._dbl) - mat(2,1) = -xy - mat(0,2) = 4._dbl*yz - mat(1,2) = -xy - mat(2,2) = 0.5_dbl*(5._dbl*zz-2._dbl*yy-1._dbl) - mat(0,3) = sqrt(2.5_dbl)*(xx-yy) - mat(1,3) = sqrt(10._dbl)*xz - mat(2,3) = -sqrt(10._dbl)*yz - mat(0,4) = sqrt(10._dbl)*xy - mat(1,4) = sqrt(10._dbl)*yz - mat(2,4) = sqrt(10._dbl)*xz - mat(0,5) = zero - mat(1,5) = sqrt(3.75_dbl)*(xx-yy) - mat(2,5) = -sqrt(15._dbl)*xy - mat(0,6) = zero - mat(1,6) = sqrt(15._dbl)*xy - mat(2,6) = sqrt(3.75_dbl)*(xx-yy) - ELSE IF (l1==7 .AND. l2==5 .AND. num==1) THEN -! d1(f,d) - mat(0,0) = sqrt(4.5_dbl)*z - mat(0,1) = -sqrt(1.5_dbl)*x - mat(0,2) = -sqrt(1.5_dbl)*y - mat(0,3) = zero - mat(0,4) = zero - mat(1,0) = sqrt(3._dbl)*x - mat(1,1) = 2._dbl*z - mat(1,2) = zero - mat(1,3) = -0.5_dbl*x - mat(1,4) = -0.5_dbl*y - mat(2,0) = sqrt(3._dbl)*y - mat(2,1) = zero - mat(2,2) = 2._dbl*z - mat(2,3) = 0.5_dbl*y - mat(2,4) = -0.5_dbl*x - mat(3,0) = zero - mat(3,1) = sqrt(2.5_dbl)*x - mat(3,2) = -sqrt(2.5_dbl)*y - mat(3,3) = sqrt(2.5_dbl)*z - mat(3,4) = zero - mat(4,0) = zero - mat(4,1) = sqrt(2.5_dbl)*y - mat(4,2) = sqrt(2.5_dbl)*x - mat(4,3) = zero - mat(4,4) = sqrt(2.5_dbl)*z - mat(5,0) = zero - mat(5,1) = zero - mat(5,2) = zero - mat(5,3) = sqrt(3.75_dbl)*x - mat(5,4) = -sqrt(3.75_dbl)*y - mat(6,0) = zero - mat(6,1) = zero - mat(6,2) = zero - mat(6,3) = sqrt(3.75_dbl)*y - mat(6,4) = sqrt(3.75_dbl)*x - ELSE IF (l1==5 .AND. l2==7 .AND. num==1) THEN -! d1(d,f) - mat(0,0) = sqrt(4.5_dbl)*z - mat(1,0) = -sqrt(1.5_dbl)*x - mat(2,0) = -sqrt(1.5_dbl)*y - mat(3,0) = zero - mat(4,0) = zero - mat(0,1) = sqrt(3._dbl)*x - mat(1,1) = 2._dbl*z - mat(2,1) = zero - mat(3,1) = -0.5_dbl*x - mat(4,1) = -0.5_dbl*y - mat(0,2) = sqrt(3._dbl)*y - mat(1,2) = zero - mat(2,2) = 2._dbl*z - mat(3,2) = 0.5_dbl*y - mat(4,2) = -0.5_dbl*x - mat(0,3) = zero - mat(1,3) = sqrt(2.5_dbl)*x - mat(2,3) = -sqrt(2.5_dbl)*y - mat(3,3) = sqrt(2.5_dbl)*z - mat(4,3) = zero - mat(0,4) = zero - mat(1,4) = sqrt(2.5_dbl)*y - mat(2,4) = sqrt(2.5_dbl)*x - mat(3,4) = zero - mat(4,4) = sqrt(2.5_dbl)*z - mat(0,5) = zero - mat(1,5) = zero - mat(2,5) = zero - mat(3,5) = sqrt(3.75_dbl)*x - mat(4,5) = -sqrt(3.75_dbl)*y - mat(0,6) = zero - mat(1,6) = zero - mat(2,6) = zero - mat(3,6) = sqrt(3.75_dbl)*y - mat(4,6) = sqrt(3.75_dbl)*x - ELSE IF (l1==7 .AND. l2==5 .AND. num==3) THEN -! d3(f,d) - mat(0,0) = sqrt(0.5_dbl)*(4._dbl*zz-3._dbl)*z - mat(0,1) = sqrt(1.5_dbl)*x*zz - mat(0,2) = sqrt(1.5_dbl)*y*zz - mat(0,3) = -sqrt(6._dbl)*(xx-yy)*z - mat(0,4) = -2._dbl*sqrt(6._dbl)*x*yz - mat(1,0) = sqrt(0.75_dbl)*(3._dbl*zz-1._dbl)*x - mat(1,1) = 0.5_dbl*(xx-5._dbl*yy)*z - mat(1,2) = 3._dbl*x*yz - mat(1,3) = (2.5_dbl*zz-xx+yy)*x - mat(1,4) = 0.5_dbl*(5._dbl*zz-4._dbl*xx)*y - mat(2,0) = sqrt(0.75_dbl)*(3._dbl*zz-1._dbl)*y - mat(2,1) = 3._dbl*x*yz - mat(2,2) = 0.5_dbl*(yy-5._dbl*xx)*z - mat(2,3) = -(2.5_dbl*zz+xx-yy)*y - mat(2,4) = 0.5_dbl*(5._dbl*zz-4._dbl*yy)*x - mat(3,0) = zero - mat(3,1) = sqrt(2.5_dbl)*(zz-2._dbl*yy)*x - mat(3,2) = -sqrt(2.5_dbl)*(zz-2._dbl*xx)*y - mat(3,3) = sqrt(2.5_dbl)*(1._dbl-2._dbl*zz)*z - mat(3,4) = zero - mat(4,0) = zero - mat(4,1) = sqrt(2.5_dbl)*(1._dbl-2._dbl*yy)*y - mat(4,2) = sqrt(2.5_dbl)*(1._dbl-2._dbl*xx)*x - mat(4,3) = zero - mat(4,4) = sqrt(2.5_dbl)*(1._dbl-2._dbl*zz)*z - mat(5,0) = sqrt(1.25_dbl)*(3._dbl*yy-xx)*x - mat(5,1) = sqrt(3.75_dbl)*(xx-yy)*z - mat(5,2) = -sqrt(15._dbl)*x*yz - mat(5,3) = -sqrt(3.75_dbl)*x*zz - mat(5,4) = sqrt(3.75_dbl)*y*zz - mat(6,0) = -sqrt(1.25_dbl)*(3._dbl*xx-yy)*y - mat(6,1) = sqrt(15._dbl)*x*yz - mat(6,2) = sqrt(3.75_dbl)*(xx-yy)*z - mat(6,3) = -sqrt(3.75_dbl)*y*zz - mat(6,4) = -sqrt(3.75_dbl)*x*zz - ELSE IF (l1==5 .AND. l2==7 .AND. num==3) THEN -! d3(d,f) - mat(0,0) = sqrt(0.5_dbl)*(4._dbl*zz-3._dbl)*z - mat(1,0) = sqrt(1.5_dbl)*x*zz - mat(2,0) = sqrt(1.5_dbl)*y*zz - mat(3,0) = -sqrt(6._dbl)*(xx-yy)*z - mat(4,0) = -2._dbl*sqrt(6._dbl)*x*yz - mat(0,1) = sqrt(0.75_dbl)*(3._dbl*zz-1._dbl)*x - mat(1,1) = 0.5_dbl*(xx-5._dbl*yy)*z - mat(2,1) = 3._dbl*x*yz - mat(3,1) = (2.5_dbl*zz-xx+yy)*x - mat(4,1) = 0.5_dbl*(5._dbl*zz-4._dbl*xx)*y - mat(0,2) = sqrt(0.75_dbl)*(3._dbl*zz-1._dbl)*y - mat(1,2) = 3._dbl*x*yz - mat(2,2) = 0.5_dbl*(yy-5._dbl*xx)*z - mat(3,2) = -(2.5_dbl*zz+xx-yy)*y - mat(4,2) = 0.5_dbl*(5._dbl*zz-4._dbl*yy)*x - mat(0,3) = zero - mat(1,3) = sqrt(2.5_dbl)*(zz-2._dbl*yy)*x - mat(2,3) = -sqrt(2.5_dbl)*(zz-2._dbl*xx)*y - mat(3,3) = sqrt(2.5_dbl)*(1._dbl-2._dbl*zz)*z - mat(4,3) = zero - mat(0,4) = zero - mat(1,4) = sqrt(2.5_dbl)*(1._dbl-2._dbl*yy)*y - mat(2,4) = sqrt(2.5_dbl)*(1._dbl-2._dbl*xx)*x - mat(3,4) = zero - mat(4,4) = sqrt(2.5_dbl)*(1._dbl-2._dbl*zz)*z - mat(0,5) = sqrt(1.25_dbl)*(3._dbl*yy-xx)*x - mat(1,5) = sqrt(3.75_dbl)*(xx-yy)*z - mat(2,5) = -sqrt(15._dbl)*x*yz - mat(3,5) = -sqrt(3.75_dbl)*x*zz - mat(4,5) = sqrt(3.75_dbl)*y*zz - mat(0,6) = -sqrt(1.25_dbl)*(3._dbl*xx-yy)*y - mat(1,6) = sqrt(15._dbl)*x*yz - mat(2,6) = sqrt(3.75_dbl)*(xx-yy)*z - mat(3,6) = -sqrt(3.75_dbl)*y*zz - mat(4,6) = -sqrt(3.75_dbl)*x*zz - ELSE IF (l1==7 .AND. l2==7 .AND. num==2) THEN -! d2(f,f) - mat(0,0) = 0.4_dbl*(3._dbl*zz-1._dbl) - mat(0,1) = sqrt(0.24_dbl)*xz - mat(0,2) = sqrt(0.24_dbl)*yz - mat(0,3) = -sqrt(0.6_dbl)*(xx-yy) - mat(0,4) = -sqrt(2.4_dbl)*xy - mat(0,5) = zero - mat(0,6) = zero - mat(1,0) = mat(0,1) - mat(2,0) = mat(0,2) - mat(3,0) = mat(0,3) - mat(4,0) = mat(0,4) - mat(5,0) = mat(0,5) - mat(6,0) = mat(0,6) - mat(1,1) = 0.3_dbl*(1._dbl-4._dbl*yy+zz) - mat(1,2) = 1.2_dbl*xy - mat(1,3) = sqrt(0.9_dbl)*xz - mat(1,4) = sqrt(0.9_dbl)*yz - mat(1,5) = -sqrt(0.15_dbl)*(xx-yy) - mat(1,6) = -sqrt(0.6_dbl)*xy - mat(2,1) = mat(1,2) - mat(3,1) = mat(1,3) - mat(4,1) = mat(1,4) - mat(5,1) = mat(1,5) - mat(6,1) = mat(1,6) - mat(2,2) = 0.3_dbl*(1._dbl-4._dbl*xx+zz) - mat(2,3) = -sqrt(0.9_dbl)*yz - mat(2,4) = sqrt(0.9_dbl)*xz - mat(2,5) = sqrt(0.6_dbl)*xy - mat(2,6) = -sqrt(0.15_dbl)*(xx-yy) - mat(3,2) = mat(2,3) - mat(4,2) = mat(2,4) - mat(5,2) = mat(2,5) - mat(6,2) = mat(2,6) - mat(3,3) = zero - mat(3,4) = zero - mat(3,5) = sqrt(1.5_dbl)*xz - mat(3,6) = sqrt(1.5_dbl)*yz - mat(4,3) = mat(3,4) - mat(5,3) = mat(3,5) - mat(6,3) = mat(3,6) - mat(4,4) = zero - mat(4,5) = -sqrt(1.5_dbl)*yz - mat(4,6) = sqrt(1.5_dbl)*xz - mat(5,4) = mat(4,5) - mat(6,4) = mat(4,6) - mat(5,5) = -0.5_dbl*(3._dbl*zz-1._dbl) - mat(5,6) = zero - mat(6,5) = zero - mat(6,6) = -0.5_dbl*(3._dbl*zz-1._dbl) - ELSE IF (l1==7 .AND. l2==7 .AND. num==4) THEN -! d4(f,f) - mat(0,0) = 0.6_dbl*(5._dbl*zz-4._dbl)*zz - mat(0,1) = sqrt(0.24_dbl)*(5._dbl*zz-2._dbl)*xz - mat(0,2) = sqrt(0.24_dbl)*(5._dbl*zz-2._dbl)*yz - mat(0,3) = -sqrt(0.6_dbl)*(xx-yy)*zz - mat(0,4) = -sqrt(2.4_dbl)*xy*zz - mat(0,5) = sqrt(3.6_dbl)*(3._dbl*yy-xx)*xz - mat(0,6) = sqrt(3.6_dbl)*(yy-3._dbl*xx)*yz - mat(1,0) = mat(0,1) - mat(2,0) = mat(0,2) - mat(3,0) = mat(0,3) - mat(4,0) = mat(0,4) - mat(5,0) = mat(0,5) - mat(6,0) = mat(0,6) - mat(1,1) = -0.4_dbl*xx - 0.5_dbl*(5._dbl*yy-3._dbl*xx)*zz - mat(1,2) = 0.4_dbl*(10._dbl*zz-1._dbl)*xy - mat(1,3) = sqrt(0.1_dbl)*(6._dbl*zz-8._dbl*yy-1._dbl)*xz - mat(1,4) = sqrt(0.1_dbl)*(2._dbl*zz-8._dbl*yy+3._dbl)*yz - mat(1,5) = sqrt(3.75_dbl)*(xx-yy-1.4_dbl*xx*xx+1.2_dbl*xx*yy+yy*yy & - ) - mat(1,6) = sqrt(0.6_dbl)*(4._dbl*zz-4._dbl*xx+1._dbl)*xy - mat(2,1) = mat(1,2) - mat(3,1) = mat(1,3) - mat(4,1) = mat(1,4) - mat(5,1) = mat(1,5) - mat(6,1) = mat(1,6) - mat(2,2) = -0.4_dbl*yy - 0.5_dbl*(5._dbl*xx-3._dbl*yy)*zz - mat(2,3) = -sqrt(0.1_dbl)*(6._dbl*zz-8._dbl*xx-1._dbl)*yz - mat(2,4) = sqrt(0.1_dbl)*(2._dbl*zz-8._dbl*xx+3._dbl)*xz - mat(2,5) = sqrt(0.6_dbl)*(4._dbl*yy-4._dbl*zz-1._dbl)*xy - mat(2,6) = sqrt(3.75_dbl)*(xx-yy-xx*xx-1.2_dbl*xx*yy+1.4_dbl*yy*yy & - ) - mat(3,2) = mat(2,3) - mat(4,2) = mat(2,4) - mat(5,2) = mat(2,5) - mat(6,2) = mat(2,6) - mat(3,3) = (2._dbl-3._dbl*zz)*zz - 4._dbl*xx*yy - mat(3,4) = 2._dbl*(xx-yy)*xy - mat(3,5) = -sqrt(1.5_dbl)*(2._dbl*zz-1._dbl)*xz - mat(3,6) = -sqrt(1.5_dbl)*(2._dbl*zz-1._dbl)*yz - mat(4,3) = mat(3,4) - mat(5,3) = mat(3,5) - mat(6,3) = mat(3,6) - mat(4,4) = -(2._dbl*zz-1._dbl)**2 + 4._dbl*xx*yy - mat(4,5) = sqrt(1.5_dbl)*(2._dbl*zz-1._dbl)*yz - mat(4,6) = -sqrt(1.5_dbl)*(2._dbl*zz-1._dbl)*xz - mat(5,4) = mat(4,5) - mat(6,4) = mat(4,6) - mat(5,5) = 1.5_dbl*(zz-1._dbl)*zz - mat(5,6) = zero - mat(6,5) = zero - mat(6,6) = 1.5_dbl*(zz-1._dbl)*zz - END IF - END IF - END SUBROUTINE set_mat -!------------------------------------------------------------------------------! - SUBROUTINE set_dmat(dmat,r,num) - IMPLICIT NONE - REAL (dbl), INTENT (IN) :: r(3) - INTEGER, INTENT (IN) :: num - REAL (dbl), INTENT (OUT) :: dmat(0:,0:,1:) - INTEGER :: l1, l2 - REAL (dbl), PARAMETER :: s3 = 0.577350269189626_dbl ! sqrt(1._dbl/3._dbl) - REAL (dbl) :: x, y, z, xx, yy, zz, xy, xz, yz - REAL (dbl) :: dc(0:6,0:6) - - l1 = size(dmat(:,0,1)) - l2 = size(dmat(0,:,1)) - x = r(1) - y = r(2) - z = r(3) - IF (l1==5 .AND. l2==3 .AND. num==1) THEN -! d1(d,p) - dmat(0:4,0:2,1:3) = zero - dmat(1,0,1) = one - dmat(0,1,1) = -s3 - dmat(3,1,1) = one - dmat(4,2,1) = one - dmat(2,0,2) = one - dmat(0,2,2) = -s3 - dmat(3,2,2) = -one - dmat(4,1,2) = one - dmat(0,0,3) = 2._dbl*s3 - dmat(1,1,3) = one - dmat(2,2,3) = one - ELSE IF (l1==3 .AND. l2==5 .AND. num==1) THEN -! d1(p,d) - dmat(0:2,0:4,1:3) = zero - dmat(0,1,1) = one - dmat(1,0,1) = -s3 - dmat(1,3,1) = one - dmat(2,4,1) = one - dmat(0,2,2) = one - dmat(2,0,2) = -s3 - dmat(2,3,2) = -one - dmat(1,4,2) = one - dmat(0,0,3) = 2._dbl*s3 - dmat(1,1,3) = one - dmat(2,2,3) = one - ELSEIF (l1==5 .AND. l2==5 .AND. num==2) THEN -! d2(d,d) - dmat(0:4,0:4,1:3) = zero - dmat(0,1,1) = s3*z - dmat(0,3,1) = -2._dbl*s3*x - dmat(0,4,1) = -2._dbl*s3*y - dmat(1,0,1) = dmat(0,1,1) - dmat(3,0,1) = dmat(0,3,1) - dmat(4,0,1) = dmat(0,4,1) - dmat(1,2,1) = y - dmat(1,3,1) = z - dmat(2,1,1) = y - dmat(3,1,1) = z - dmat(2,2,1) = -2._dbl*x - dmat(2,4,1) = z - dmat(4,2,1) = z - dmat(0,2,2) = s3*z - dmat(0,3,2) = 2._dbl*s3*y - dmat(0,4,2) = -2._dbl*s3*x - dmat(2,0,2) = dmat(0,2,2) - dmat(3,0,2) = dmat(0,3,2) - dmat(4,0,2) = dmat(0,4,2) - dmat(1,1,2) = -2._dbl*y - dmat(1,2,2) = x - dmat(1,4,2) = z - dmat(2,1,2) = x - dmat(4,1,2) = z - dmat(2,3,2) = -z - dmat(3,2,2) = -z - dmat(0,0,3) = 2._dbl*z - dmat(0,1,3) = s3*x - dmat(0,2,3) = s3*y - dmat(1,0,3) = dmat(0,1,3) - dmat(2,0,3) = dmat(0,2,3) - dmat(1,3,3) = x - dmat(1,4,3) = y - dmat(3,1,3) = x - dmat(4,1,3) = y - dmat(2,3,3) = -y - dmat(2,4,3) = x - dmat(3,2,3) = -y - dmat(4,2,3) = x - dmat(3,3,3) = -2._dbl*z - dmat(4,4,3) = -2._dbl*z - ELSE IF (l1==7 .AND. l2==3 .AND. num==2) THEN -! d2(f,p) - dmat(0,0,1) = zero - dmat(0,1,1) = -sqrt(6._dbl)*z - dmat(0,2,1) = zero - dmat(1,0,1) = 4._dbl*z - dmat(1,1,1) = -2._dbl*x - dmat(1,2,1) = -y - dmat(2,0,1) = zero - dmat(2,1,1) = -y - dmat(2,2,1) = zero - dmat(3,0,1) = 2._dbl*sqrt(2.5_dbl)*x - dmat(3,1,1) = sqrt(10._dbl)*z - dmat(3,2,1) = zero - dmat(4,0,1) = sqrt(10._dbl)*y - dmat(4,1,1) = zero - dmat(4,2,1) = sqrt(10._dbl)*z - dmat(5,0,1) = zero - dmat(5,1,1) = sqrt(15._dbl)*x - dmat(5,2,1) = -sqrt(15._dbl)*y - dmat(6,0,1) = zero - dmat(6,1,1) = sqrt(15._dbl)*y - dmat(6,2,1) = sqrt(15._dbl)*x - - dmat(0,0,2) = zero - dmat(0,1,2) = zero - dmat(0,2,2) = -sqrt(6._dbl)*z - dmat(1,0,2) = zero - dmat(1,1,2) = zero - dmat(1,2,2) = -x - dmat(2,0,2) = 4._dbl*z - dmat(2,1,2) = -x - dmat(2,2,2) = -2._dbl*y - dmat(3,0,2) = -2._dbl*sqrt(2.5_dbl)*y - dmat(3,1,2) = zero - dmat(3,2,2) = -sqrt(10._dbl)*z - dmat(4,0,2) = sqrt(10._dbl)*x - dmat(4,1,2) = sqrt(10._dbl)*z - dmat(4,2,2) = zero - dmat(5,0,2) = zero - dmat(5,1,2) = -sqrt(15._dbl)*y - dmat(5,2,2) = -sqrt(15._dbl)*x - dmat(6,0,2) = zero - dmat(6,1,2) = sqrt(15._dbl)*x - dmat(6,2,2) = -sqrt(15._dbl)*y - - dmat(0,0,3) = 6._dbl*sqrt(1.5_dbl)*z - dmat(0,1,3) = -sqrt(6._dbl)*x - dmat(0,2,3) = -sqrt(6._dbl)*y - dmat(1,0,3) = 4._dbl*x - dmat(1,1,3) = 5._dbl*z - dmat(1,2,3) = zero - dmat(2,0,3) = 4._dbl*y - dmat(2,1,3) = zero - dmat(2,2,3) = 5._dbl*z - dmat(3,0,3) = zero - dmat(3,1,3) = sqrt(10._dbl)*x - dmat(3,2,3) = -sqrt(10._dbl)*y - dmat(4,0,3) = zero - dmat(4,1,3) = sqrt(10._dbl)*y - dmat(4,2,3) = sqrt(10._dbl)*x - dmat(5,0,3) = zero - dmat(5,1,3) = zero - dmat(5,2,3) = zero - dmat(6,0,3) = zero - dmat(6,1,3) = zero - dmat(6,2,3) = zero - ELSE IF (l1==3 .AND. l2==7 .AND. num==2) THEN -! d2(p,f) - dmat(0,0,1) = zero - dmat(1,0,1) = -sqrt(6._dbl)*z - dmat(2,0,1) = zero - dmat(0,1,1) = 4._dbl*z - dmat(1,1,1) = -2._dbl*x - dmat(2,1,1) = -y - dmat(0,2,1) = zero - dmat(1,2,1) = -y - dmat(2,2,1) = zero - dmat(0,3,1) = 2._dbl*sqrt(2.5_dbl)*x - dmat(1,3,1) = sqrt(10._dbl)*z - dmat(2,3,1) = zero - dmat(0,4,1) = sqrt(10._dbl)*y - dmat(1,4,1) = zero - dmat(2,4,1) = sqrt(10._dbl)*z - dmat(0,5,1) = zero - dmat(1,5,1) = sqrt(15._dbl)*x - dmat(2,5,1) = -sqrt(15._dbl)*y - dmat(0,6,1) = zero - dmat(1,6,1) = sqrt(15._dbl)*y - dmat(2,6,1) = sqrt(15._dbl)*x - - dmat(0,0,2) = zero - dmat(1,0,2) = zero - dmat(2,0,2) = -sqrt(6._dbl)*z - dmat(0,1,2) = zero - dmat(1,1,2) = zero - dmat(2,1,2) = -x - dmat(0,2,2) = 4._dbl*z - dmat(1,2,2) = -x - dmat(2,2,2) = -2._dbl*y - dmat(0,3,2) = -2._dbl*sqrt(2.5_dbl)*y - dmat(1,3,2) = zero - dmat(2,3,2) = -sqrt(10._dbl)*z - dmat(0,4,2) = sqrt(10._dbl)*x - dmat(1,4,2) = sqrt(10._dbl)*z - dmat(2,4,2) = zero - dmat(0,5,2) = zero - dmat(1,5,2) = -sqrt(15._dbl)*y - dmat(2,5,2) = -sqrt(15._dbl)*x - dmat(0,6,2) = zero - dmat(1,6,2) = sqrt(15._dbl)*x - dmat(2,6,2) = -sqrt(15._dbl)*y - - dmat(0,0,3) = 6._dbl*sqrt(1.5_dbl)*z - dmat(1,0,3) = -sqrt(6._dbl)*x - dmat(2,0,3) = -sqrt(6._dbl)*y - dmat(0,1,3) = 4._dbl*x - dmat(1,1,3) = 5._dbl*z - dmat(2,1,3) = zero - dmat(0,2,3) = 4._dbl*y - dmat(1,2,3) = zero - dmat(2,2,3) = 5._dbl*z - dmat(0,3,3) = zero - dmat(1,3,3) = sqrt(10._dbl)*x - dmat(2,3,3) = -sqrt(10._dbl)*y - dmat(0,4,3) = zero - dmat(1,4,3) = sqrt(10._dbl)*y - dmat(2,4,3) = sqrt(10._dbl)*x - dmat(0,5,3) = zero - dmat(1,5,3) = zero - dmat(2,5,3) = zero - dmat(0,6,3) = zero - dmat(1,6,3) = zero - dmat(2,6,3) = zero - ELSE IF (l1==7 .AND. l2==5 .AND. num==1) THEN -! d1(f,d) - dmat(0:6,0:4,1:3) = zero - dmat(0,1,1) = -sqrt(1.5_dbl) - dmat(1,0,1) = sqrt(3._dbl) - dmat(1,3,1) = -0.5_dbl - dmat(2,4,1) = -0.5_dbl - dmat(3,1,1) = sqrt(2.5_dbl) - dmat(4,2,1) = sqrt(2.5_dbl) - dmat(5,3,1) = sqrt(3.75_dbl) - dmat(6,4,1) = sqrt(3.75_dbl) - - dmat(0,2,2) = -sqrt(1.5_dbl) - dmat(1,4,2) = -0.5_dbl - dmat(2,0,2) = sqrt(3._dbl) - dmat(2,3,2) = 0.5_dbl - dmat(3,2,2) = -sqrt(2.5_dbl) - dmat(4,1,2) = sqrt(2.5_dbl) - dmat(5,4,2) = -sqrt(3.75_dbl) - dmat(6,3,2) = sqrt(3.75_dbl) - - dmat(0,0,3) = sqrt(4.5_dbl) - dmat(1,1,3) = 2._dbl - dmat(2,2,3) = 2._dbl - dmat(3,3,3) = sqrt(2.5_dbl) - dmat(4,4,3) = sqrt(2.5_dbl) - ELSE IF (l1==5 .AND. l2==7 .AND. num==1) THEN -! d1(d,f) - dmat(0:4,0:6,1:3) = zero - dmat(1,0,1) = -sqrt(1.5_dbl) - dmat(0,1,1) = sqrt(3._dbl) - dmat(3,1,1) = -0.5_dbl - dmat(4,2,1) = -0.5_dbl - dmat(1,3,1) = sqrt(2.5_dbl) - dmat(2,4,1) = sqrt(2.5_dbl) - dmat(3,5,1) = sqrt(3.75_dbl) - dmat(4,6,1) = sqrt(3.75_dbl) - - dmat(2,0,2) = -sqrt(1.5_dbl) - dmat(4,1,2) = -0.5_dbl - dmat(0,2,2) = sqrt(3._dbl) - dmat(3,2,2) = 0.5_dbl - dmat(2,3,2) = -sqrt(2.5_dbl) - dmat(1,4,2) = sqrt(2.5_dbl) - dmat(4,5,2) = -sqrt(3.75_dbl) - dmat(3,6,2) = sqrt(3.75_dbl) - - dmat(0,0,3) = sqrt(4.5_dbl) - dmat(1,1,3) = 2._dbl - dmat(2,2,3) = 2._dbl - dmat(3,3,3) = sqrt(2.5_dbl) - dmat(4,4,3) = sqrt(2.5_dbl) - ELSE IF (l1==7 .AND. l2==5 .AND. num==3) THEN - xx = x*x - yy = y*y - zz = z*z - xy = x*y - xz = x*z - yz = y*z -! d3(f,d) - dmat(0,0,1) = zero - dmat(0,1,1) = sqrt(1.5_dbl)*zz - dmat(0,2,1) = zero - dmat(0,3,1) = -2._dbl*sqrt(6._dbl)*xz - dmat(0,4,1) = -2._dbl*sqrt(6._dbl)*yz - dmat(1,0,1) = sqrt(0.75_dbl)*(3._dbl*zz-1._dbl) - dmat(1,1,1) = xz - dmat(1,2,1) = 3._dbl*yz - dmat(1,3,1) = 2.5_dbl*zz-3._dbl*xx+yy - dmat(1,4,1) = -4._dbl*xy - dmat(2,0,1) = zero - dmat(2,1,1) = 3._dbl*yz - dmat(2,2,1) = -5._dbl*xz - dmat(2,3,1) = -2._dbl*xy - dmat(2,4,1) = 0.5_dbl*(5._dbl*zz-4._dbl*yy) - dmat(3,0,1) = zero - dmat(3,1,1) = sqrt(2.5_dbl)*(zz-2._dbl*yy) - dmat(3,2,1) = sqrt(2.5_dbl)*4._dbl*xy - dmat(3,3,1) = zero - dmat(3,4,1) = zero - dmat(4,0,1) = zero - dmat(4,1,1) = zero - dmat(4,2,1) = sqrt(2.5_dbl)*(1._dbl-6._dbl*xx) - dmat(4,3,1) = zero - dmat(4,4,1) = zero - dmat(5,0,1) = 3._dbl*sqrt(1.25_dbl)*(yy-xx) - dmat(5,1,1) = sqrt(15._dbl)*xz - dmat(5,2,1) = -sqrt(15._dbl)*yz - dmat(5,3,1) = -sqrt(3.75_dbl)*zz - dmat(5,4,1) = zero - dmat(6,0,1) = -sqrt(5._dbl)*3._dbl*xy - dmat(6,1,1) = sqrt(15._dbl)*yz - dmat(6,2,1) = sqrt(15._dbl)*xz - dmat(6,3,1) = zero - dmat(6,4,1) = -sqrt(3.75_dbl)*zz - - dmat(0,0,2) = zero - dmat(0,1,2) = zero - dmat(0,2,2) = sqrt(1.5_dbl)*zz - dmat(0,3,2) = 2._dbl*sqrt(6._dbl)*yz - dmat(0,4,2) = -2._dbl*sqrt(6._dbl)*xz - dmat(1,0,2) = zero - dmat(1,1,2) = -5._dbl*yz - dmat(1,2,2) = 3._dbl*xz - dmat(1,3,2) = 2._dbl*xy - dmat(1,4,2) = 0.5_dbl*(5._dbl*zz-4._dbl*xx) - dmat(2,0,2) = sqrt(0.75_dbl)*(3._dbl*zz-1._dbl) - dmat(2,1,2) = 3._dbl*xz - dmat(2,2,2) = yz - dmat(2,3,2) = -(2.5_dbl*zz+xx-3._dbl*yy) - dmat(2,4,2) = -4._dbl*xy - dmat(3,0,2) = zero - dmat(3,1,2) = -sqrt(2.5_dbl)*4._dbl*xy - dmat(3,2,2) = -sqrt(2.5_dbl)*(zz-2._dbl*xx) - dmat(3,3,2) = zero - dmat(3,4,2) = zero - dmat(4,0,2) = zero - dmat(4,1,2) = sqrt(2.5_dbl)*(1._dbl-6._dbl*yy) - dmat(4,2,2) = zero - dmat(4,3,2) = zero - dmat(4,4,2) = zero - dmat(5,0,2) = 3._dbl*sqrt(5._dbl)*xy - dmat(5,1,2) = -sqrt(15._dbl)*yz - dmat(5,2,2) = -sqrt(15._dbl)*xz - dmat(5,3,2) = zero - dmat(5,4,2) = sqrt(3.75_dbl)*zz - dmat(6,0,2) = -sqrt(1.25_dbl)*(3._dbl*xx-3._dbl*yy) - dmat(6,1,2) = sqrt(15._dbl)*xz - dmat(6,2,2) = -sqrt(15._dbl)*yz - dmat(6,3,2) = -sqrt(3.75_dbl)*zz - dmat(6,4,2) = zero - - dmat(0,0,3) = sqrt(0.5_dbl)*(12._dbl*zz-3._dbl) - dmat(0,1,3) = 2._dbl*sqrt(1.5_dbl)*xz - dmat(0,2,3) = 2._dbl*sqrt(1.5_dbl)*yz - dmat(0,3,3) = -sqrt(6._dbl)*(xx-yy) - dmat(0,4,3) = -2._dbl*sqrt(6._dbl)*xy - dmat(1,0,3) = 3._dbl*sqrt(3._dbl)*xz - dmat(1,1,3) = 0.5_dbl*(xx-5._dbl*yy) - dmat(1,2,3) = 3._dbl*xy - dmat(1,3,3) = 5._dbl*xz - dmat(1,4,3) = 5._dbl*yz - dmat(2,0,3) = 3._dbl*sqrt(3._dbl)*yz - dmat(2,1,3) = 3._dbl*xy - dmat(2,2,3) = 0.5_dbl*(yy-5._dbl*xx) - dmat(2,3,3) = -5._dbl*yz - dmat(2,4,3) = 5._dbl*xz - dmat(3,0,3) = zero - dmat(3,1,3) = 2._dbl*sqrt(2.5_dbl)*xz - dmat(3,2,3) = -2._dbl*sqrt(2.5_dbl)*yz - dmat(3,3,3) = sqrt(2.5_dbl)*(1._dbl-6._dbl*zz) - dmat(3,4,3) = zero - dmat(4,0,3) = zero - dmat(4,1,3) = zero - dmat(4,2,3) = zero - dmat(4,3,3) = zero - dmat(4,4,3) = sqrt(2.5_dbl)*(1._dbl-6._dbl*zz) - dmat(5,0,3) = zero - dmat(5,1,3) = sqrt(3.75_dbl)*(xx-yy) - dmat(5,2,3) = -sqrt(15._dbl)*xy - dmat(5,3,3) = -sqrt(15._dbl)*xz - dmat(5,4,3) = sqrt(15._dbl)*yz - dmat(6,0,3) = zero - dmat(6,1,3) = sqrt(15._dbl)*xy - dmat(6,2,3) = sqrt(3.75_dbl)*(xx-yy) - dmat(6,3,3) = -sqrt(15._dbl)*yz - dmat(6,4,3) = -sqrt(15._dbl)*xz - - ELSE IF (l1==5 .AND. l2==7 .AND. num==3) THEN - xx = x*x - yy = y*y - zz = z*z - xy = x*y - xz = x*z - yz = y*z -! d3(d,f) - dmat(0,0,1) = zero - dmat(1,0,1) = sqrt(1.5_dbl)*zz - dmat(2,0,1) = zero - dmat(3,0,1) = -2._dbl*sqrt(6._dbl)*xz - dmat(4,0,1) = -2._dbl*sqrt(6._dbl)*yz - dmat(0,1,1) = sqrt(0.75_dbl)*(3._dbl*zz-1._dbl) - dmat(1,1,1) = xz - dmat(2,1,1) = 3._dbl*yz - dmat(3,1,1) = 2.5_dbl*zz-3._dbl*xx+yy - dmat(4,1,1) = -4._dbl*xy - dmat(0,2,1) = zero - dmat(1,2,1) = 3._dbl*yz - dmat(2,2,1) = -5._dbl*xz - dmat(3,2,1) = -2._dbl*xy - dmat(4,2,1) = 0.5_dbl*(5._dbl*zz-4._dbl*yy) - dmat(0,3,1) = zero - dmat(1,3,1) = sqrt(2.5_dbl)*(zz-2._dbl*yy) - dmat(2,3,1) = sqrt(2.5_dbl)*4._dbl*xy - dmat(3,3,1) = zero - dmat(4,3,1) = zero - dmat(0,4,1) = zero - dmat(1,4,1) = zero - dmat(2,4,1) = sqrt(2.5_dbl)*(1._dbl-6._dbl*xx) - dmat(3,4,1) = zero - dmat(4,4,1) = zero - dmat(0,5,1) = 3._dbl*sqrt(1.25_dbl)*(yy-xx) - dmat(1,5,1) = sqrt(15._dbl)*xz - dmat(2,5,1) = -sqrt(15._dbl)*yz - dmat(3,5,1) = -sqrt(3.75_dbl)*zz - dmat(4,5,1) = zero - dmat(0,6,1) = -sqrt(5._dbl)*3._dbl*xy - dmat(1,6,1) = sqrt(15._dbl)*yz - dmat(2,6,1) = sqrt(15._dbl)*xz - dmat(3,6,1) = zero - dmat(4,6,1) = -sqrt(3.75_dbl)*zz - - dmat(0,0,2) = zero - dmat(1,0,2) = zero - dmat(2,0,2) = sqrt(1.5_dbl)*zz - dmat(3,0,2) = 2._dbl*sqrt(6._dbl)*yz - dmat(4,0,2) = -2._dbl*sqrt(6._dbl)*xz - dmat(0,1,2) = zero - dmat(1,1,2) = -5._dbl*yz - dmat(2,1,2) = 3._dbl*xz - dmat(3,1,2) = 2._dbl*xy - dmat(4,1,2) = 0.5_dbl*(5._dbl*zz-4._dbl*xx) - dmat(0,2,2) = sqrt(0.75_dbl)*(3._dbl*zz-1._dbl) - dmat(1,2,2) = 3._dbl*xz - dmat(2,2,2) = yz - dmat(3,2,2) = -(2.5_dbl*zz+xx-3._dbl*yy) - dmat(4,2,2) = -4._dbl*xy - dmat(0,3,2) = zero - dmat(1,3,2) = -sqrt(2.5_dbl)*4._dbl*xy - dmat(2,3,2) = -sqrt(2.5_dbl)*(zz-2._dbl*xx) - dmat(3,3,2) = zero - dmat(4,3,2) = zero - dmat(0,4,2) = zero - dmat(1,4,2) = sqrt(2.5_dbl)*(1._dbl-6._dbl*yy) - dmat(2,4,2) = zero - dmat(3,4,2) = zero - dmat(4,4,2) = zero - dmat(0,5,2) = sqrt(5._dbl)*3._dbl*xy - dmat(1,5,2) = -sqrt(15._dbl)*yz - dmat(2,5,2) = -sqrt(15._dbl)*xz - dmat(3,5,2) = zero - dmat(4,5,2) = sqrt(3.75_dbl)*zz - dmat(0,6,2) = -sqrt(1.25_dbl)*(3._dbl*xx-3._dbl*yy) - dmat(1,6,2) = sqrt(15._dbl)*xz - dmat(2,6,2) = -sqrt(15._dbl)*yz - dmat(3,6,2) = -sqrt(3.75_dbl)*zz - dmat(4,6,2) = zero - - dmat(0,0,3) = sqrt(0.5_dbl)*(12._dbl*zz-3._dbl) - dmat(1,0,3) = 2._dbl*sqrt(1.5_dbl)*xz - dmat(2,0,3) = 2._dbl*sqrt(1.5_dbl)*yz - dmat(3,0,3) = -sqrt(6._dbl)*(xx-yy) - dmat(4,0,3) = -2._dbl*sqrt(6._dbl)*xy - dmat(0,1,3) = 3._dbl*sqrt(3._dbl)*xz - dmat(1,1,3) = 0.5_dbl*(xx-5._dbl*yy) - dmat(2,1,3) = 3._dbl*xy - dmat(3,1,3) = 5._dbl*xz - dmat(4,1,3) = 5._dbl*yz - dmat(0,2,3) = 3._dbl*sqrt(3._dbl)*yz - dmat(1,2,3) = 3._dbl*xy - dmat(2,2,3) = 0.5_dbl*(yy-5._dbl*xx) - dmat(3,2,3) = -5._dbl*yz - dmat(4,2,3) = 5._dbl*xz - dmat(0,3,3) = zero - dmat(1,3,3) = 2._dbl*sqrt(2.5_dbl)*xz - dmat(2,3,3) = -2._dbl*sqrt(2.5_dbl)*yz - dmat(3,3,3) = sqrt(2.5_dbl)*(1._dbl-6._dbl*zz) - dmat(4,3,3) = zero - dmat(0,4,3) = zero - dmat(1,4,3) = zero - dmat(2,4,3) = zero - dmat(3,4,3) = zero - dmat(4,4,3) = sqrt(2.5_dbl)*(1._dbl-6._dbl*zz) - dmat(0,5,3) = zero - dmat(1,5,3) = sqrt(3.75_dbl)*(xx-yy) - dmat(2,5,3) = -sqrt(15._dbl)*xy - dmat(3,5,3) = -sqrt(15._dbl)*xz - dmat(4,5,3) = sqrt(15._dbl)*yz - dmat(0,6,3) = zero - dmat(1,6,3) = sqrt(15._dbl)*xy - dmat(2,6,3) = sqrt(3.75_dbl)*(xx-yy) - dmat(3,6,3) = -sqrt(15._dbl)*yz - dmat(4,6,3) = -sqrt(15._dbl)*xz - ELSE IF (l1==7 .AND. l2==7 .AND. num==2) THEN -! d2(f,f) - dmat(0,0,1) = zero - dmat(0,1,1) = sqrt(0.24_dbl)*z - dmat(0,2,1) = zero - dmat(0,3,1) = -2._dbl*sqrt(0.6_dbl)*x - dmat(0,4,1) = -sqrt(2.4_dbl)*y - dmat(0,5,1) = zero - dmat(0,6,1) = zero - dmat(1,0,1) = dmat(0,1,1) - dmat(2,0,1) = dmat(0,2,1) - dmat(3,0,1) = dmat(0,3,1) - dmat(4,0,1) = dmat(0,4,1) - dmat(5,0,1) = dmat(0,5,1) - dmat(6,0,1) = dmat(0,6,1) - dmat(1,1,1) = zero - dmat(1,2,1) = 1.2_dbl*y - dmat(1,3,1) = sqrt(0.9_dbl)*z - dmat(1,4,1) = zero - dmat(1,5,1) = -2._dbl*sqrt(0.15_dbl)*x - dmat(1,6,1) = -sqrt(0.6_dbl)*y - dmat(2,1,1) = dmat(1,2,1) - dmat(3,1,1) = dmat(1,3,1) - dmat(4,1,1) = dmat(1,4,1) - dmat(5,1,1) = dmat(1,5,1) - dmat(6,1,1) = dmat(1,6,1) - dmat(2,2,1) = -2.4_dbl*x - dmat(2,3,1) = zero - dmat(2,4,1) = sqrt(0.9_dbl)*z - dmat(2,5,1) = sqrt(0.6_dbl)*y - dmat(2,6,1) = -2._dbl*sqrt(0.15_dbl)*x - dmat(3,2,1) = dmat(2,3,1) - dmat(4,2,1) = dmat(2,4,1) - dmat(5,2,1) = dmat(2,5,1) - dmat(6,2,1) = dmat(2,6,1) - dmat(3,3,1) = zero - dmat(3,4,1) = zero - dmat(3,5,1) = sqrt(1.5_dbl)*z - dmat(3,6,1) = zero - dmat(4,3,1) = dmat(3,4,1) - dmat(5,3,1) = dmat(3,5,1) - dmat(6,3,1) = dmat(3,6,1) - dmat(4,4,1) = zero - dmat(4,5,1) = zero - dmat(4,6,1) = sqrt(1.5_dbl)*z - dmat(5,4,1) = dmat(4,5,1) - dmat(6,4,1) = dmat(4,6,1) - dmat(5,5,1) = zero - dmat(5,6,1) = zero - dmat(6,5,1) = zero - dmat(6,6,1) = zero - - dmat(0,0,2) = zero - dmat(0,1,2) = zero - dmat(0,2,2) = sqrt(0.24_dbl)*z - dmat(0,3,2) = 2._dbl*sqrt(0.6_dbl)*y - dmat(0,4,2) = -sqrt(2.4_dbl)*x - dmat(0,5,2) = zero - dmat(0,6,2) = zero - dmat(1,0,2) = dmat(0,1,2) - dmat(2,0,2) = dmat(0,2,2) - dmat(3,0,2) = dmat(0,3,2) - dmat(4,0,2) = dmat(0,4,2) - dmat(5,0,2) = dmat(0,5,2) - dmat(6,0,2) = dmat(0,6,2) - dmat(1,1,2) = -2.4_dbl*y - dmat(1,2,2) = 1.2_dbl*x - dmat(1,3,2) = zero - dmat(1,4,2) = sqrt(0.9_dbl)*z - dmat(1,5,2) = 2._dbl*sqrt(0.15_dbl)*y - dmat(1,6,2) = -sqrt(0.6_dbl)*x - dmat(2,1,2) = dmat(1,2,2) - dmat(3,1,2) = dmat(1,3,2) - dmat(4,1,2) = dmat(1,4,2) - dmat(5,1,2) = dmat(1,5,2) - dmat(6,1,2) = dmat(1,6,2) - dmat(2,2,2) = zero - dmat(2,3,2) = -sqrt(0.9_dbl)*z - dmat(2,4,2) = zero - dmat(2,5,2) = sqrt(0.6_dbl)*x - dmat(2,6,2) = 2._dbl*sqrt(0.15_dbl)*y - dmat(3,2,2) = dmat(2,3,2) - dmat(4,2,2) = dmat(2,4,2) - dmat(5,2,2) = dmat(2,5,2) - dmat(6,2,2) = dmat(2,6,2) - dmat(3,3,2) = zero - dmat(3,4,2) = zero - dmat(3,5,2) = zero - dmat(3,6,2) = sqrt(1.5_dbl)*z - dmat(4,3,2) = dmat(3,4,2) - dmat(5,3,2) = dmat(3,5,2) - dmat(6,3,2) = dmat(3,6,2) - dmat(4,4,2) = zero - dmat(4,5,2) = -sqrt(1.5_dbl)*z - dmat(4,6,2) = zero - dmat(5,4,2) = dmat(4,5,2) - dmat(6,4,2) = dmat(4,6,2) - dmat(5,5,2) = zero - dmat(5,6,2) = zero - dmat(6,5,2) = zero - dmat(6,6,2) = zero - - dmat(0,0,3) = 2.4_dbl*z - dmat(0,1,3) = sqrt(0.24_dbl)*x - dmat(0,2,3) = sqrt(0.24_dbl)*y - dmat(0,3,3) = zero - dmat(0,4,3) = zero - dmat(0,5,3) = zero - dmat(0,6,3) = zero - dmat(1,0,3) = dmat(0,1,3) - dmat(2,0,3) = dmat(0,2,3) - dmat(3,0,3) = dmat(0,3,3) - dmat(4,0,3) = dmat(0,4,3) - dmat(5,0,3) = dmat(0,5,3) - dmat(6,0,3) = dmat(0,6,3) - dmat(1,1,3) = 0.6_dbl*z - dmat(1,2,3) = zero - dmat(1,3,3) = sqrt(0.9_dbl)*x - dmat(1,4,3) = sqrt(0.9_dbl)*y - dmat(1,5,3) = zero - dmat(1,6,3) = zero - dmat(2,1,3) = dmat(1,2,3) - dmat(3,1,3) = dmat(1,3,3) - dmat(4,1,3) = dmat(1,4,3) - dmat(5,1,3) = dmat(1,5,3) - dmat(6,1,3) = dmat(1,6,3) - dmat(2,2,3) = 0.6_dbl*z - dmat(2,3,3) = -sqrt(0.9_dbl)*y - dmat(2,4,3) = sqrt(0.9_dbl)*x - dmat(2,5,3) = zero - dmat(2,6,3) = zero - dmat(3,2,3) = dmat(2,3,3) - dmat(4,2,3) = dmat(2,4,3) - dmat(5,2,3) = dmat(2,5,3) - dmat(6,2,3) = dmat(2,6,3) - dmat(3,3,3) = zero - dmat(3,4,3) = zero - dmat(3,5,3) = sqrt(1.5_dbl)*x - dmat(3,6,3) = sqrt(1.5_dbl)*y - dmat(4,3,3) = dmat(3,4,3) - dmat(5,3,3) = dmat(3,5,3) - dmat(6,3,3) = dmat(3,6,3) - dmat(4,4,3) = zero - dmat(4,5,3) = -sqrt(1.5_dbl)*y - dmat(4,6,3) = sqrt(1.5_dbl)*x - dmat(5,4,3) = dmat(4,5,3) - dmat(6,4,3) = dmat(4,6,3) - dmat(5,5,3) = -3._dbl*z - dmat(5,6,3) = zero - dmat(6,5,3) = zero - dmat(6,6,3) = -3._dbl*z - ELSE IF (l1==7 .AND. l2==7 .AND. num==4) THEN - xx = x*x - yy = y*y - zz = z*z - xy = x*y - xz = x*z - yz = y*z -! d4(f,f) - dmat(0,0,1) = zero - dmat(0,1,1) = sqrt(0.24_dbl)*(5._dbl*zz-2._dbl)*z - dmat(0,2,1) = zero - dmat(0,3,1) = -sqrt(0.6_dbl)*2._dbl*x*zz - dmat(0,4,1) = -sqrt(2.4_dbl)*y*zz - dmat(0,5,1) = sqrt(3.6_dbl)*(3._dbl*yy-3._dbl*xx)*z - dmat(0,6,1) = -sqrt(3.6_dbl)*6._dbl*x*yz - dmat(1,0,1) = dmat(0,1,1) - dmat(2,0,1) = dmat(0,2,1) - dmat(3,0,1) = dmat(0,3,1) - dmat(4,0,1) = dmat(0,4,1) - dmat(5,0,1) = dmat(0,5,1) - dmat(6,0,1) = dmat(0,6,1) - dmat(1,1,1) = -0.8_dbl*x + 3._dbl*x*zz - dmat(1,2,1) = 0.4_dbl*(10._dbl*zz-1._dbl)*y - dmat(1,3,1) = sqrt(0.1_dbl)*(6._dbl*zz-8._dbl*yy-1._dbl)*z - dmat(1,4,1) = zero - dmat(1,5,1) = sqrt(3.75_dbl)*(2._dbl*x-5.6_dbl*x*xx+2.4_dbl*x*yy) - dmat(1,6,1) = sqrt(0.6_dbl)*(4._dbl*zz-12._dbl*xx+1._dbl)*y - dmat(2,1,1) = dmat(1,2,1) - dmat(3,1,1) = dmat(1,3,1) - dmat(4,1,1) = dmat(1,4,1) - dmat(5,1,1) = dmat(1,5,1) - dmat(6,1,1) = dmat(1,6,1) - dmat(2,2,1) = -5._dbl*x*zz - dmat(2,3,1) = sqrt(0.1_dbl)*16._dbl*x*yz - dmat(2,4,1) = sqrt(0.1_dbl)*(2._dbl*zz-24._dbl*xx+3._dbl)*z - dmat(2,5,1) = sqrt(0.6_dbl)*(4._dbl*yy-4._dbl*zz-1._dbl)*y - dmat(2,6,1) = sqrt(3.75_dbl)*(2._dbl*x-4._dbl*x*xx-2.4_dbl*x*yy) - dmat(3,2,1) = dmat(2,3,1) - dmat(4,2,1) = dmat(2,4,1) - dmat(5,2,1) = dmat(2,5,1) - dmat(6,2,1) = dmat(2,6,1) - dmat(3,3,1) = -8._dbl*x*yy - dmat(3,4,1) = 2._dbl*(3._dbl*xx-yy)*y - dmat(3,5,1) = -sqrt(1.5_dbl)*(2._dbl*zz-1._dbl)*z - dmat(3,6,1) = zero - dmat(4,3,1) = dmat(3,4,1) - dmat(5,3,1) = dmat(3,5,1) - dmat(6,3,1) = dmat(3,6,1) - dmat(4,4,1) = 8._dbl*x*yy - dmat(4,5,1) = zero - dmat(4,6,1) = -sqrt(1.5_dbl)*(2._dbl*zz-1._dbl)*z - dmat(5,4,1) = dmat(4,5,1) - dmat(6,4,1) = dmat(4,6,1) - dmat(5,5,1) = zero - dmat(5,6,1) = zero - dmat(6,5,1) = zero - dmat(6,6,1) = zero - - dmat(0,0,2) = zero - dmat(0,1,2) = zero - dmat(0,2,2) = sqrt(0.24_dbl)*(5._dbl*zz-2._dbl)*z - dmat(0,3,2) = sqrt(0.6_dbl)*2._dbl*y*zz - dmat(0,4,2) = -sqrt(2.4_dbl)*x*zz - dmat(0,5,2) = sqrt(3.6_dbl)*6._dbl*y*xz - dmat(0,6,2) = sqrt(3.6_dbl)*(3._dbl*yy-3._dbl*xx)*z - dmat(1,0,2) = dmat(0,1,2) - dmat(2,0,2) = dmat(0,2,2) - dmat(3,0,2) = dmat(0,3,2) - dmat(4,0,2) = dmat(0,4,2) - dmat(5,0,2) = dmat(0,5,2) - dmat(6,0,2) = dmat(0,6,2) - dmat(1,1,2) = -5._dbl*y*zz - dmat(1,2,2) = 0.4_dbl*(10._dbl*zz-1._dbl)*x - dmat(1,3,2) = sqrt(0.1_dbl)*(-16._dbl*y)*xz - dmat(1,4,2) = sqrt(0.1_dbl)*(2._dbl*zz-24._dbl*yy+3._dbl)*z - dmat(1,5,2) = sqrt(3.75_dbl)*(-2._dbl*y+2.4_dbl*xx*y+4._dbl*y*yy) - dmat(1,6,2) = sqrt(0.6_dbl)*(4._dbl*zz-4._dbl*xx+1._dbl)*x - dmat(2,1,2) = dmat(1,2,2) - dmat(3,1,2) = dmat(1,3,2) - dmat(4,1,2) = dmat(1,4,2) - dmat(5,1,2) = dmat(1,5,2) - dmat(6,1,2) = dmat(1,6,2) - dmat(2,2,2) = -0.8_dbl*y + 3._dbl*y*zz - dmat(2,3,2) = -sqrt(0.1_dbl)*(6._dbl*zz-8._dbl*xx-1._dbl)*z - dmat(2,4,2) = zero - dmat(2,5,2) = sqrt(0.6_dbl)*(12._dbl*yy-4._dbl*zz-1._dbl)*x - dmat(2,6,2) = sqrt(3.75_dbl)*(-2._dbl*y-2.4_dbl*xx*y+5.6_dbl*y*yy) - dmat(3,2,2) = dmat(2,3,2) - dmat(4,2,2) = dmat(2,4,2) - dmat(5,2,2) = dmat(2,5,2) - dmat(6,2,2) = dmat(2,6,2) - dmat(3,3,2) = -8._dbl*xx*y - dmat(3,4,2) = 2._dbl*(xx-3._dbl*yy)*x - dmat(3,5,2) = zero - dmat(3,6,2) = -sqrt(1.5_dbl)*(2._dbl*zz-1._dbl)*z - dmat(4,3,2) = dmat(3,4,2) - dmat(5,3,2) = dmat(3,5,2) - dmat(6,3,2) = dmat(3,6,2) - dmat(4,4,2) = 8._dbl*xx*y - dmat(4,5,2) = sqrt(1.5_dbl)*(2._dbl*zz-1._dbl)*z - dmat(4,6,2) = zero - dmat(5,4,2) = dmat(4,5,2) - dmat(6,4,2) = dmat(4,6,2) - dmat(5,5,2) = zero - dmat(5,6,2) = zero - dmat(6,5,2) = zero - dmat(6,6,2) = zero - - dmat(0,0,3) = 0.6_dbl*(20._dbl*zz-8._dbl)*z - dmat(0,1,3) = sqrt(0.24_dbl)*(15._dbl*zz-2._dbl)*x - dmat(0,2,3) = sqrt(0.24_dbl)*(15._dbl*zz-2._dbl)*y - dmat(0,3,3) = -sqrt(0.6_dbl)*(xx-yy)*2._dbl*z - dmat(0,4,3) = -sqrt(2.4_dbl)*xy*2._dbl*z - dmat(0,5,3) = sqrt(3.6_dbl)*(3._dbl*yy-xx)*x - dmat(0,6,3) = sqrt(3.6_dbl)*(yy-3._dbl*xx)*y - dmat(1,0,3) = dmat(0,1,3) - dmat(2,0,3) = dmat(0,2,3) - dmat(3,0,3) = dmat(0,3,3) - dmat(4,0,3) = dmat(0,4,3) - dmat(5,0,3) = dmat(0,5,3) - dmat(6,0,3) = dmat(0,6,3) - dmat(1,1,3) = -(5._dbl*yy-3._dbl*xx)*z - dmat(1,2,3) = 8._dbl*z*xy - dmat(1,3,3) = sqrt(0.1_dbl)*(18._dbl*zz-8._dbl*yy-1._dbl)*x - dmat(1,4,3) = sqrt(0.1_dbl)*(6._dbl*zz-8._dbl*yy+3._dbl)*y - dmat(1,5,3) = zero - dmat(1,6,3) = sqrt(0.6_dbl)*8._dbl*z*xy - dmat(2,1,3) = dmat(1,2,3) - dmat(3,1,3) = dmat(1,3,3) - dmat(4,1,3) = dmat(1,4,3) - dmat(5,1,3) = dmat(1,5,3) - dmat(6,1,3) = dmat(1,6,3) - dmat(2,2,3) = -(5._dbl*xx-3._dbl*yy)*z - dmat(2,3,3) = -sqrt(0.1_dbl)*(18._dbl*zz-8._dbl*xx-1._dbl)*y - dmat(2,4,3) = sqrt(0.1_dbl)*(6._dbl*zz-8._dbl*xx+3._dbl)*x - dmat(2,5,3) = sqrt(0.6_dbl)*(-8._dbl*z)*xy - dmat(2,6,3) = zero - dmat(3,2,3) = dmat(2,3,3) - dmat(4,2,3) = dmat(2,4,3) - dmat(5,2,3) = dmat(2,5,3) - dmat(6,2,3) = dmat(2,6,3) - dmat(3,3,3) = (4._dbl-12._dbl*zz)*z - dmat(3,4,3) = zero - dmat(3,5,3) = -sqrt(1.5_dbl)*(6._dbl*zz-1._dbl)*x - dmat(3,6,3) = -sqrt(1.5_dbl)*(6._dbl*zz-1._dbl)*y - dmat(4,3,3) = dmat(3,4,3) - dmat(5,3,3) = dmat(3,5,3) - dmat(6,3,3) = dmat(3,6,3) - dmat(4,4,3) = -2._dbl*(2._dbl*zz-1._dbl)*4._dbl*z - dmat(4,5,3) = sqrt(1.5_dbl)*(6._dbl*zz-1._dbl)*y - dmat(4,6,3) = -sqrt(1.5_dbl)*(6._dbl*zz-1._dbl)*x - dmat(5,4,3) = dmat(4,5,3) - dmat(6,4,3) = dmat(4,6,3) - dmat(5,5,3) = 1.5_dbl*(4._dbl*zz-2._dbl)*z - dmat(5,6,3) = zero - dmat(6,5,3) = zero - dmat(6,6,3) = 1.5_dbl*(4._dbl*zz-2._dbl)*z - END IF - dc(0:l1-1,0:l2-1) = r(1)*dmat(:,:,1) + r(2)*dmat(:,:,2) + & - r(3)*dmat(:,:,3) - dmat(:,:,1) = dmat(:,:,1) - r(1)*dc(0:l1-1,0:l2-1) - dmat(:,:,2) = dmat(:,:,2) - r(2)*dc(0:l1-1,0:l2-1) - dmat(:,:,3) = dmat(:,:,3) - r(3)*dc(0:l1-1,0:l2-1) - END SUBROUTINE set_dmat -!------------------------------------------------------------------------------! - END MODULE slater_koster_matr -!------------------------------------------------------------------------------! diff --git a/src/slater_koster_matr.c b/src/slater_koster_matr.c new file mode 100644 index 0000000..31dd696 --- /dev/null +++ b/src/slater_koster_matr.c @@ -0,0 +1,1458 @@ +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright (C) 2000 CP2K developers group */ +/*----------------------------------------------------------------------------*/ +//> Calculation of two-center s-f Slater-Koster integrals +// A.K. McMahan, Phys. Rev. B, 58, p4293 (1998) +// +// Sign convention: t(lp,lq)=(-1)**(lp+lq)*(lp m,lq m) +// where (lp m,lq m) is the integral in the diatomic coordinate system +// +//< + +//#include +#include + +#include "slater_koster.h" + +void gmat(int l1, int l2, double r[3], double gmu[7][7][4], double dll[7][7][7]) { + int mu, mo, m1, m2, i; + double gsl[7][4]; + + mu = (l1 < l2) ? l1 : l2; + mo = (l1 > l2) ? l1 : l2; + m1 = 2 * l1; + m2 = 2 * l2; + + memset(gsl, 0, sizeof(gsl)); + for (i = 0; i <= mo; i++) { + sph(i, r, gsl[i]); + } + + memset(dll, 0, sizeof(dll)); + for (i = 0; i <= 2 * mo; i++) { + dll[i][i][0] = 1.0; + } + + if (l1 == 0) { + // s-l + for (i = 0; i <= m2; i++) { + gmu[0][i][0] = gsl[i][l2]; + } + } else if (l2 == 0) { + // l-s + for (i = 0; i <= m1; i++) { + gmu[i][0][0] = gsl[i][l1]; + } + } else { + switch (l1) { + case 1: + switch (l2) { + case 1: + // p-p + out_prod(&dll[0][2][2], &gsl[1], &gsl[1]); + memcpy(gmu[0][0], &dll[2][0][0], sizeof(gmu[0][0])); + memcpy(gmu[0][1], (-dll[2][0][2] + dll[2][0][0]), sizeof(gmu[0][1])); + break; + case 2: + // p-d + out_prod(&dll[0][4][3], gsl[1], gsl[2]); + set_mat(dll[0][4][1], r, 1); + memcpy(gmu[0][0], &dll[3][0][0], sizeof(gmu[0][0])); + memcpy(gmu[0][1], -2.0 / sqrt(3.0) * &dll[3][0][3] + dll[3][0][1], sizeof(gmu[0][1])); + break; + case 3: + // p-f + out_prod(&dll[0][6][4], gsl[1], gsl[3]); + set_mat(dll[0][6][2], r, 2); + memcpy(gmu[0][0], dll[4][0][0], sizeof(gmu[0][0])); + memcpy(gmu[0][1], -sqrt(1.5) * dll[4][0][4] + 0.5 * dll[4][0][2], sizeof(gmu[0][1])); + break; + } + break; + case 2: + switch (l2) { + case 1: + // d-p + out_prod(dll[0][2][3], gsl[2], gsl[1]); + set_mat(dll[0][2][1], r, 1); + memcpy(gmu[0][0], dll[3][0][0], sizeof(gmu[0][0])); + memcpy(gmu[0][1], -2.0 / sqrt(3.0) * dll[3][0][3] + dll[3][0][1], sizeof(gmu[0][1])); + break; + case 2: + // d-d + out_prod(dll[0][4][4], gsl[2], gsl[2]); + set_mat(dll[0][4][2], r, 2); + memcpy(gmu[0][0], dll[4][0][0], sizeof(gmu[0][0])); + memcpy(gmu[0][1], -4.0 / 3.0 * dll[4][0][4] + dll[4][0][2] + dll[4][0][0], sizeof(gmu[0][1])); + memcpy(gmu[0][2], 1.0 / 3.0 * dll[4][0][4] - dll[4][0][2], sizeof(gmu[0][2])); + break; + case 3: + // d-f + out_prod(dll[0][6][5], gsl[2], gsl[3]); + set_mat(dll[0][6][3], r, 3); + set_mat(dll[0][6][1], r, 1); + memcpy(gmu[0][0], dll[5][0][0], sizeof(gmu[0][0])); + memcpy(gmu[0][1], -sqrt(2.0) * dll[5][0][5] + 0.5 * dll[5][0][3] + 0.5 * dll[5][0][1], sizeof(gmu[0][1])); + memcpy(gmu[0][2], sqrt(0.2) * dll[5][0][5] - sqrt(0.4) * dll[5][0][3], sizeof(gmu[0][2])); + break; + } + break; + case 3: + switch (l2) { + case 1: + // f-p + out_prod(dll[0][2][4], gsl[3], gsl[1]); + set_mat(dll[0][2][2], r, 2); + memcpy(gmu[0][0], dll[4][0][0], sizeof(gmu[0][0])); + memcpy(gmu[0][1], -sqrt(1.5) * dll[4][0][4] + 0.5 * dll[4][0][2], sizeof(gmu[0][1])); + break; + case 2: + // f-d + out_prod(dll[0][4][5], gsl[3], gsl[2]); + set_mat(dll[0][4][3], r, 3); + set_mat(dll[0][4][1], r, 1); + memcpy(gmu[0][0], dll[5][0][0], sizeof(gmu[0][0])); + memcpy(gmu[0][1], -sqrt(2.0) * dll[5][0][5] + 0.5 * dll[5][0][3] + 0.5 * dll[5][0][1], sizeof(gmu[0][1])); + memcpy(gmu[0][2], sqrt(0.2) * dll[5][0][5] - sqrt(0.4) * dll[5][0][3], sizeof(gmu[0][2])); + break; + case 3: + // f-f + out_prod(dll[0][6][6], gsl[3], gsl[3]); + set_mat(dll[0][6][4], r, 4); + set_mat(dll[0][6][2], r, 2); + memcpy(gmu[0][0], dll[6][0][0], sizeof(gmu[0][0])); + memcpy(gmu[0][1], -1.5 * dll[6][0][6] + 0.625 * dll[6][0][4] + 0.625 * dll[6][0][2] + 0.625 * dll[6][0][0], sizeof(gmu[0][1])); + memcpy(gmu[0][2], 0.6 * dll[6][0][6] - dll[6][0][4], sizeof(gmu[0][2])); + memcpy(gmu[0][3], -0.1 * dll[6][0][6] + 0.375 * dll[6][0][4] - 0.625 * dll[6][0][2] + 0.375 * dll[6][0][0], sizeof(gmu[0][3])); + break; + } + break; + } + } +} + +/*----------------------------------------------------------------------------*/ +void dgmat(int l1, int l2, double r[3], double dgmu[7][7][4][3], double ddll[7][7][7][3]) { + int mu, mo, m1, m2, i; + double gsl[7][4]; + double dgsl[7][4][3]; + + mu = (l1 < l2) ? l1 : l2; + mo = (l1 > l2) ? l1 : l2; + m1 = 2 * l1; + m2 = 2 * l2; + memset(gsl, 0, sizeof(gsl)); + memset(dgsl, 0, sizeof(dgsl)); + + for (i = 0; i <= mo; i++) { + sph(i, r, gsl[i]); + dsph(i, r, dgsl[i]); + } + + if (l1 == 0) { + // s-l + memcpy(dgmu[0][0], dgsl[l2], sizeof(dgmu[0][0])); + } else if (l2 == 0) { + // l-s + memcpy(dgmu[0][0], dgsl[l1], sizeof(dgmu[0][0])); + } else { + switch (l1) { + case 1: + switch (l2) { + case 1: + // p-p + out_dprod(ddll[0][0][2][0], gsl[1], gsl[1], dgsl[1][0], dgsl[1][0]); + out_dprod(ddll[0][0][2][1], gsl[1], gsl[1], dgsl[1][1], dgsl[1][1]); + out_dprod(ddll[0][0][2][2], gsl[1], gsl[1], dgsl[1][2], dgsl[1][2]); + memcpy(dgmu[0][0][0], ddll[0][0][2], sizeof(dgmu[0][0][0])); + memcpy(dgmu[0][0][1], ddll[0][0][2], sizeof(dgmu[0][0][1])); + break; + case 2: + // p-d + out_dprod(&ddll[0][0][3][0], gsl[1], gsl[2], dgsl[1][0], dgsl[2][0]); + out_dprod(&ddll[0][0][3][1], gsl[1], gsl[2], dgsl[1][1], dgsl[2][1]); + out_dprod(&ddll[0][0][3][2], gsl[1], gsl[2], dgsl[1][2], dgsl[2][2]); + set_dmat(ddll[0][0][1], r, 1); + memcpy(dgmu[0][0][0], ddll[0][0][3], sizeof(dgmu[0][0][0])); + memcpy(dgmu[0][0][1], ddll[0][0][3], sizeof(dgmu[0][0][1])); + break; + case 3: + // p-f + out_dprod(&ddll[0][0][4][0], gsl[1], gsl[3], dgsl[1][0], dgsl[3][0]); + out_dprod(&ddll[0][0][4][1], gsl[1], gsl[3], dgsl[1][1], dgsl[3][1]); + out_dprod(&ddll[0][0][4][2], gsl[1], gsl[3], dgsl[1][2], dgsl[3][2]); + set_dmat(ddll[0][0][2], r, 2); + memcpy(dgmu[0][0][0], ddll[0][0][4], sizeof(dgmu[0][0][0])); + memcpy(dgmu[0][0][1], ddll[0][0][4], sizeof(dgmu[0][0][1])); + break; + } + break; + case 2: + switch (l2) { + case 1: + // d-p + out_dprod(&ddll[0][0][3][0], &gsl[2], gsl[1], dgsl[2][0], dgsl[1][0]); + out_dprod(&ddll[0][0][3][1], gsl[2], gsl[1], dgsl[2][1], dgsl[1][1]); + out_dprod(&ddll[0][0][3][2], gsl[2], gsl[1], dgsl[2][2], dgsl[1][2]); + set_dmat(ddll[0][0][1], r, 1); + memcpy(dgmu[0][0][0], ddll[0][0][3], sizeof(dgmu[0][0][0])); + memcpy(dgmu[0][0][1], ddll[0][0][3], sizeof(dgmu[0][0][1])); + break; + case 2: + // d-d + out_dprod(&ddll[0][0][4][0], gsl[2], gsl[2], dgsl[2][0], dgsl[2][0]); + out_dprod(&ddll[0][0][4][1], gsl[2], gsl[2], dgsl[2][1], dgsl[2][1]); + out_dprod(&ddll[0][0][4][2], gsl[2], gsl[2], dgsl[2][2], dgsl[2][2]); + set_dmat(ddll[0][0][2], r, 2); + memcpy(dgmu[0][0][0], ddll[0][0][4], sizeof(dgmu[0][0][0])); + memcpy(dgmu[0][0][1], ddll[0][0][4], sizeof(dgmu[0][0][1])); + memcpy(dgmu[0][0][2], ddll[0][0][4], sizeof(dgmu[0][0][2])); + break; + case 3: + // d-f + out_dprod(&ddll[0][0][5][0], gsl[2], gsl[3], dgsl[2][0], dgsl[3][0]); + out_dprod(&ddll[0][0][5][1], gsl[2], gsl[3], dgsl[2][1], dgsl[3][1]); + out_dprod(&ddll[0][0][5][2], gsl[2], gsl[3], dgsl[2][2], dgsl[3][2]); + set_dmat(ddll[0][0][3], r, 3); + set_dmat(ddll[0][0][1], r, 1); + memcpy(dgmu[0][0][0], ddll[0][0][5], sizeof(dgmu[0][0][0])); + memcpy(dgmu[0][0][1], ddll[0][0][5], sizeof(dgmu[0][0][1])); + memcpy(dgmu[0][0][2], ddll[0][0][5], sizeof(dgmu[0][0][2])); + break; + } + break; + case 3: + switch (l2) { + case 1: + // f-p + out_dprod(&ddll[0][0][4][0], gsl[3], gsl[1], dgsl[3][0], dgsl[1][0]); + out_dprod(&ddll[0][0][4][1], gsl[3], gsl[1], dgsl[3][1], dgsl[1][1]); + out_dprod(&ddll[0][0][4][2], gsl[3], gsl[1], dgsl[3][2], dgsl[1][2]); + set_dmat(ddll[0][0][2], r, 2); + memcpy(dgmu[0][0][0], ddll[0][0][4], sizeof(dgmu[0][0][0])); + memcpy(dgmu[0][0][1], ddll[0][0][4], sizeof(dgmu[0][0][1])); + break; + case 2: + // f-d + out_dprod(&ddll[0][0][5][0], gsl[3], gsl[2], dgsl[3][0], dgsl[2][0]); + out_dprod(&ddll[0][0][5][1], gsl[3], gsl[2], dgsl[3][1], dgsl[2][1]); + out_dprod(&ddll[0][0][5][2], gsl[3], gsl[2], dgsl[3][2], dgsl[2][2]); + set_dmat(ddll[0][0][3], r, 3); + set_dmat(ddll[0][0][1], r, 1); + memcpy(dgmu[0][0][0], ddll[0][0][5], sizeof(dgmu[0][0][0])); + memcpy(dgmu[0][0][1], ddll[0][0][5], sizeof(dgmu[0][0][1])); + memcpy(dgmu[0][0][2], ddll[0][0][5], sizeof(dgmu[0][0][2])); + break; + case 3: + // f-f + out_dprod(&ddll[0][0][6][0], gsl[3], gsl[3], dgsl[3][0], dgsl[3][0]); + out_dprod(&ddll[0][0][6][1], gsl[3], gsl[3], dgsl[3][1], dgsl[3][1]); + out_dprod(&ddll[0][0][6][2], gsl[3], gsl[3], dgsl[3][2], dgsl[3][2]); + set_dmat(ddll[0][0][4], r, 4); + set_dmat(ddll[0][0][2], r, 2); + memcpy(dgmu[0][0][0], ddll[0][0][6], sizeof(dgmu[0][0][0])); + memcpy(dgmu[0][0][1], ddll[0][0][6], sizeof(dgmu[0][0][1])); + memcpy(dgmu[0][0][2], ddll[0][0][6], sizeof(dgmu[0][0][2])); + memcpy(dgmu[0][0][3], ddll[0][0][6], sizeof(dgmu[0][0][3])); + break; + } + break; + } + } +} +/*----------------------------------------------------------------------------*/ +void set_mat(double mat[][3], double r[3], int num) { + + int l1, l2; + const double s3 = 0.577350269189626; // sqrt(1.0/3.0) + double x, y, z, xx, yy, zz, xy, xz, yz; + + l1 = sizeof(mat) / sizeof(mat[0]); + l2 = sizeof(mat[0]) / sizeof(mat[0][0]); + x = r[0]; + y = r[1]; + z = r[2]; + + if (l1 == 5 && l2 == 3 && num == 1) { +// d1(d,p) + mat[0][0] = 2.0 * s3 * z; + mat[1][0] = x; + mat[2][0] = y; + mat[0][1] = -s3 * x; + mat[1][1] = z; + mat[3][1] = x; + mat[4][1] = y; + mat[0][2] = -s3 * y; + mat[2][2] = z; + mat[3][2] = -y; + mat[4][2] = x; + } else if (l1 == 3 && l2 == 5 && num == 1) { +// d1(p,d) + mat[0][0] = 2.0 * s3 * z; + mat[0][1] = x; + mat[0][2] = y; + mat[1][0] = -s3 * x; + mat[1][1] = z; + mat[1][3] = x; + mat[1][4] = y; + mat[2][0] = -s3 * y; + mat[2][2] = z; + mat[2][3] = -y; + mat[2][4] = x; + } else { + xx = x * x; + yy = y * y; + zz = z * z; + xy = x * y; + xz = x * z; + yz = y * z; + if (l1 == 5 && l2 == 5 && num == 2) { +// d2(d,d) + mat[0][0] = zz - 2.0/3.0; + mat[0][1] = s3*xz; + mat[0][2] = s3*yz; + mat[0][3] = -s3*(xx-yy); + mat[0][4] = -2.0*s3*xy; + mat[1][0] = mat[0][1]; + mat[2][0] = mat[0][2]; + mat[3][0] = mat[0][3]; + mat[4][0] = mat[0][4]; + mat[1][1] = -yy; + mat[1][2] = xy; + mat[1][3] = xz; + mat[1][4] = yz; + mat[2][1] = mat[1][2]; + mat[3][1] = mat[1][3]; + mat[4][1] = mat[1][4]; + mat[2][2] = -xx; + mat[2][3] = -yz; + mat[2][4] = xz; + mat[3][2] = mat[2][3]; + mat[4][2] = mat[2][4]; + mat[3][3] = -zz; + mat[3][4] = 0.0; + mat[4][3] = 0.0; + mat[4][4] = -zz; + } else if (l1 == 7 && l2 == 3 && num == 2) { +// d2(f,p) + mat[0][0] = sqrt(1.50)*(3.0*zz-1.0); + mat[0][1] = -sqrt(6.0)*xz; + mat[0][2] = -sqrt(6.0)*yz; + mat[1][0] = 4.0*xz; + mat[1][1] = 0.50*(5.0*zz-2.0*xx-1.0); + mat[1][2] = -xy; + mat[2][0] = 4.0*yz; + mat[2][1] = -xy; + mat[2][2] = 0.50*(5.0*zz-2.0*yy-1.0); + mat[3][0] = sqrt(2.50)*(xx-yy); + mat[3][1] = sqrt(10.0)*xz; + mat[3][2] = -sqrt(10.0)*yz; + mat[4][0] = sqrt(10.0)*xy; + mat[4][1] = sqrt(10.0)*yz; + mat[4][2] = sqrt(10.0)*xz; + mat[5][0] = 0.0; + mat[5][1] = sqrt(3.750)*(xx-yy); + mat[5][2] = -sqrt(15.0)*xy; + mat[6][0] = 0.0; + mat[6][1] = sqrt(15.0)*xy; + mat[6][2] = sqrt(3.750)*(xx-yy); + } else if (l1 == 3 && l2 == 7 && num == 2) { +// d2(p,f) + mat[0][0] = sqrt(1.50)*(3.0*zz-1.0); + mat[1][0] = -sqrt(6.0)*xz; + mat[2][0] = -sqrt(6.0)*yz; + mat[0][1] = 4.0*xz; + mat[1][1] = 0.50*(5.0*zz-2.0*xx-1.0); + mat[2][1] = -xy; + mat[0][2] = 4.0*yz; + mat[1][2] = -xy; + mat[2][2] = 0.50*(5.0*zz-2.0*yy-1.0); + mat[0][3] = sqrt(2.50)*(xx-yy); + mat[1][3] = sqrt(10.0)*xz; + mat[2][3] = -sqrt(10.0)*yz; + mat[0][4] = sqrt(10.0)*xy; + mat[1][4] = sqrt(10.0)*yz; + mat[2][4] = sqrt(10.0)*xz; + mat[0][5] = 0.0; + mat[1][5] = sqrt(3.750)*(xx-yy); + mat[2][5] = -sqrt(15.0)*xy; + mat[0][6] = 0.0; + mat[1][6] = sqrt(15.0)*xy; + mat[2][6] = sqrt(3.750)*(xx-yy); + } else if (l1 == 7 && l2 == 5 && num == 1) { +// d1(f,d) + mat[0][0] = sqrt(4.50)*z; + mat[0][1] = -sqrt(1.50)*x; + mat[0][2] = -sqrt(1.50)*y; + mat[0][3] = 0.0; + mat[0][4] = 0.0; + mat[1][0] = sqrt(3.0)*x; + mat[1][1] = 2.0*z; + mat[1][2] = 0.0; + mat[1][3] = -0.50*x; + mat[1][4] = -0.50*y; + mat[2][0] = sqrt(3.0)*y; + mat[2][1] = 0.0; + mat[2][2] = 2.0*z; + mat[2][3] = 0.50*y; + mat[2][4] = -0.50*x; + mat[3][0] = 0.0; + mat[3][1] = sqrt(2.50)*x; + mat[3][2] = -sqrt(2.50)*y; + mat[3][3] = sqrt(2.50)*z; + mat[3][4] = 0.0; + mat[4][0] = 0.0; + mat[4][1] = sqrt(2.50)*y; + mat[4][2] = sqrt(2.50)*x; + mat[4][3] = 0.0; + mat[4][4] = sqrt(2.50)*z; + mat[5][0] = 0.0; + mat[5][1] = 0.0; + mat[5][2] = 0.0; + mat[5][3] = sqrt(3.750)*x; + mat[5][4] = -sqrt(3.750)*y; + mat[6][0] = 0.0; + mat[6][1] = 0.0; + mat[6][2] = 0.0; + mat[6][3] = sqrt(3.750)*y; + mat[6][4] = sqrt(3.750)*x; + } else if (l1 == 5 && l2 == 7 && num == 1) { +// d1(d,f) + mat[0][0] = sqrt(4.50)*z; + mat[1][0] = -sqrt(1.50)*x; + mat[2][0] = -sqrt(1.50)*y; + mat[3][0] = 0.0; + mat[4][0] = 0.0; + mat[0][1] = sqrt(3.0)*x; + mat[1][1] = 2.0*z; + mat[2][1] = 0.0; + mat[3][1] = -0.50*x; + mat[4][1] = -0.50*y; + mat[0][2] = sqrt(3.0)*y; + mat[1][2] = 0.0; + mat[2][2] = 2.0*z; + mat[3][2] = 0.50*y; + mat[4][2] = -0.50*x; + mat[0][3] = 0.0; + mat[1][3] = sqrt(2.50)*x; + mat[2][3] = -sqrt(2.50)*y; + mat[3][3] = sqrt(2.50)*z; + mat[4][3] = 0.0; + mat[0][4] = 0.0; + mat[1][4] = sqrt(2.50)*y; + mat[2][4] = sqrt(2.50)*x; + mat[3][4] = 0.0; + mat[4][4] = sqrt(2.50)*z; + mat[0][5] = 0.0; + mat[1][5] = 0.0; + mat[2][5] = 0.0; + mat[3][5] = sqrt(3.750)*x; + mat[4][5] = -sqrt(3.750)*y; + mat[0][6] = 0.0; + mat[1][6] = 0.0; + mat[2][6] = 0.0; + mat[3][6] = sqrt(3.750)*y; + mat[4][6] = sqrt(3.750)*x; + } else if (l1 == 7 && l2 == 5 && num == 3) { +// d3(f,d) + mat[0][0] = sqrt(0.50)*(4.0*zz-3.0)*z; + mat[0][1] = sqrt(1.50)*x*zz; + mat[0][2] = sqrt(1.50)*y*zz; + mat[0][3] = -sqrt(6.0)*(xx-yy)*z; + mat[0][4] = -2.0*sqrt(6.0)*x*yz; + mat[1][0] = sqrt(0.750)*(3.0*zz-1.0)*x; + mat[1][1] = 0.50*(xx-5.0*yy)*z; + mat[1][2] = 3.0*x*yz; + mat[1][3] = (2.50*zz-xx+yy)*x; + mat[1][4] = 0.50*(5.0*zz-4.0*xx)*y; + mat[2][0] = sqrt(0.750)*(3.0*zz-1.0)*y; + mat[2][1] = 3.0*x*yz; + mat[2][2] = 0.50*(yy-5.0*xx)*z; + mat[2][3] = -(2.50*zz+xx-yy)*y; + mat[2][4] = 0.50*(5.0*zz-4.0*yy)*x; + mat[3][0] = 0.0; + mat[3][1] = sqrt(2.50)*(zz-2.0*yy)*x; + mat[3][2] = -sqrt(2.50)*(zz-2.0*xx)*y; + mat[3][3] = sqrt(2.50)*(1.0-2.0*zz)*z; + mat[3][4] = 0.0; + mat[4][0] = 0.0; + mat[4][1] = sqrt(2.50)*(1.0-2.0*yy)*y; + mat[4][2] = sqrt(2.50)*(1.0-2.0*xx)*x; + mat[4][3] = 0.0; + mat[4][4] = sqrt(2.50)*(1.0-2.0*zz)*z; + mat[5][0] = sqrt(1.250)*(3.0*yy-xx)*x; + mat[5][1] = sqrt(3.750)*(xx-yy)*z; + mat[5][2] = -sqrt(15.0)*x*yz; + mat[5][3] = -sqrt(3.750)*x*zz; + mat[5][4] = sqrt(3.750)*y*zz; + mat[6][0] = -sqrt(1.250)*(3.0*xx-yy)*y; + mat[6][1] = sqrt(15.0)*x*yz; + mat[6][2] = sqrt(3.750)*(xx-yy)*z; + mat[6][3] = -sqrt(3.750)*y*zz; + mat[6][4] = -sqrt(3.750)*x*zz; + } else if (l1 == 5 && l2 == 7 && num == 3) { +// d3(d,f) + mat[0][0] = sqrt(0.50)*(4.0*zz-3.0)*z; + mat[1][0] = sqrt(1.50)*x*zz; + mat[2][0] = sqrt(1.50)*y*zz; + mat[3][0] = -sqrt(6.0)*(xx-yy)*z; + mat[4][0] = -2.0*sqrt(6.0)*x*yz; + mat[0][1] = sqrt(0.750)*(3.0*zz-1.0)*x; + mat[1][1] = 0.50*(xx-5.0*yy)*z; + mat[2][1] = 3.0*x*yz; + mat[3][1] = (2.50*zz-xx+yy)*x; + mat[4][1] = 0.50*(5.0*zz-4.0*xx)*y; + mat[0][2] = sqrt(0.750)*(3.0*zz-1.0)*y; + mat[1][2] = 3.0*x*yz; + mat[2][2] = 0.50*(yy-5.0*xx)*z; + mat[3][2] = -(2.50*zz+xx-yy)*y; + mat[4][2] = 0.50*(5.0*zz-4.0*yy)*x; + mat[0][3] = 0.0; + mat[1][3] = sqrt(2.50)*(zz-2.0*yy)*x; + mat[2][3] = -sqrt(2.50)*(zz-2.0*xx)*y; + mat[3][3] = sqrt(2.50)*(1.0-2.0*zz)*z; + mat[4][3] = 0.0; + mat[0][4] = 0.0; + mat[1][4] = sqrt(2.50)*(1.0-2.0*yy)*y; + mat[2][4] = sqrt(2.50)*(1.0-2.0*xx)*x; + mat[3][4] = 0.0; + mat[4][4] = sqrt(2.50)*(1.0-2.0*zz)*z; + mat[0][5] = sqrt(1.250)*(3.0*yy-xx)*x; + mat[1][5] = sqrt(3.750)*(xx-yy)*z; + mat[2][5] = -sqrt(15.0)*x*yz; + mat[3][5] = -sqrt(3.750)*x*zz; + mat[4][5] = sqrt(3.750)*y*zz; + mat[0][6] = -sqrt(1.250)*(3.0*xx-yy)*y; + mat[1][6] = sqrt(15.0)*x*yz; + mat[2][6] = sqrt(3.750)*(xx-yy)*z; + mat[3][6] = -sqrt(3.750)*y*zz; + mat[4][6] = -sqrt(3.750)*x*zz; + } else if (l1 == 7 && l2 == 7 && num == 2) { +// d2(f,f) + mat[0][0] = 0.40*(3.0*zz-1.0); + mat[0][1] = sqrt(0.240)*xz; + mat[0][2] = sqrt(0.240)*yz; + mat[0][3] = -sqrt(0.60)*(xx-yy); + mat[0][4] = -sqrt(2.40)*xy; + mat[0][5] = 0.0; + mat[0][6] = 0.0; + mat[1][0] = mat[0][1]; + mat[2][0] = mat[0][2]; + mat[3][0] = mat[0][3]; + mat[4][0] = mat[0][4]; + mat[5][0] = mat[0][5]; + mat[6][0] = mat[0][6]; + mat[1][1] = 0.30*(1.0-4.0*yy+zz); + mat[1][2] = 1.20*xy; + mat[1][3] = sqrt(0.90)*xz; + mat[1][4] = sqrt(0.90)*yz; + mat[1][5] = -sqrt(0.150)*(xx-yy); + mat[1][6] = -sqrt(0.60)*xy; + mat[2][1] = mat[1][2]; + mat[3][1] = mat[1][3]; + mat[4][1] = mat[1][4]; + mat[5][1] = mat[1][5]; + mat[6][1] = mat[1][6]; + mat[2][2] = 0.30*(1.0-4.0*xx+zz); + mat[2][3] = -sqrt(0.90)*yz; + mat[2][4] = sqrt(0.90)*xz; + mat[2][5] = sqrt(0.60)*xy; + mat[2][6] = -sqrt(0.150)*(xx-yy); + mat[3][2] = mat[2][3]; + mat[4][2] = mat[2][4]; + mat[5][2] = mat[2][5]; + mat[6][2] = mat[2][6]; + mat[3][3] = 0.0; + mat[3][4] = 0.0; + mat[3][5] = sqrt(1.50)*xz; + mat[3][6] = sqrt(1.50)*yz; + mat[4][3] = mat[3][4]; + mat[5][3] = mat[3][5]; + mat[6][3] = mat[3][6]; + mat[4][4] = 0.0; + mat[4][5] = -sqrt(1.50)*yz; + mat[4][6] = sqrt(1.50)*xz; + mat[5][4] = mat[4][5]; + mat[6][4] = mat[4][6]; + mat[5][5] = -0.50*(3.0*zz-1.0); + mat[5][6] = 0.0; + mat[6][5] = 0.0; + mat[6][6] = -0.50*(3.0*zz-1.0); + } else if (l1 == 7 && l2 == 7 && num == 4) { +// d4(f,f) + mat[0][0] = 0.60*(5.0*zz-4.0)*zz; + mat[0][1] = sqrt(0.240)*(5.0*zz-2.0)*xz; + mat[0][2] = sqrt(0.240)*(5.0*zz-2.0)*yz; + mat[0][3] = -sqrt(0.60)*(xx-yy)*zz; + mat[0][4] = -sqrt(2.40)*xy*zz; + mat[0][5] = sqrt(3.60)*(3.0*yy-xx)*xz; + mat[0][6] = sqrt(3.60)*(yy-3.0*xx)*yz; + mat[1][0] = mat[0][1]; + mat[2][0] = mat[0][2]; + mat[3][0] = mat[0][3]; + mat[4][0] = mat[0][4]; + mat[5][0] = mat[0][5]; + mat[6][0] = mat[0][6]; + mat[1][1] = -0.40*xx - 0.50*(5.0*yy-3.0*xx)*zz; + mat[1][2] = 0.40*(10.0*zz-1.0)*xy; + mat[1][3] = sqrt(0.10)*(6.0*zz-8.0*yy-1.0)*xz; + mat[1][4] = sqrt(0.10)*(2.0*zz-8.0*yy+3.0)*yz; + mat[1][5] = sqrt(3.750)*(xx-yy-1.40*xx*xx+1.20*xx*yy+yy*yy); + mat[1][6] = sqrt(0.60)*(4.0*zz-4.0*xx+1.0)*xy; + mat[2][1] = mat[1][2]; + mat[3][1] = mat[1][3]; + mat[4][1] = mat[1][4]; + mat[5][1] = mat[1][5]; + mat[6][1] = mat[1][6]; + mat[2][2] = -0.40*yy - 0.50*(5.0*xx-3.0*yy)*zz; + mat[2][3] = -sqrt(0.10)*(6.0*zz-8.0*xx-1.0)*yz; + mat[2][4] = sqrt(0.10)*(2.0*zz-8.0*xx+3.0)*xz; + mat[2][5] = sqrt(0.60)*(4.0*yy-4.0*zz-1.0)*xy; + mat[2][6] = sqrt(3.750)*(xx-yy-xx*xx-1.20*xx*yy+1.40*yy*yy); + mat[3][2] = mat[2][3]; + mat[4][2] = mat[2][4]; + mat[5][2] = mat[2][5]; + mat[6][2] = mat[2][6]; + mat[3][3] = (2.0-3.0*zz)*zz - 4.0*xx*yy; + mat[3][4] = 2.0*(xx-yy)*xy; + mat[3][5] = -sqrt(1.50)*(2.0*zz-1.0)*xz; + mat[3][6] = -sqrt(1.50)*(2.0*zz-1.0)*yz; + mat[4][3] = mat[3][4]; + mat[5][3] = mat[3][5]; + mat[6][3] = mat[3][6]; + mat[4][4] = -(2.0*zz-1.0)*(2.0*zz-1.0) + 4.0*xx*yy; + mat[4][5] = sqrt(1.50)*(2.0*zz-1.0)*yz; + mat[4][6] = -sqrt(1.50)*(2.0*zz-1.0)*xz; + mat[5][4] = mat[4][5]; + mat[6][4] = mat[4][6]; + mat[5][5] = 1.50*(zz-1.0)*zz; + mat[5][6] = 0.0; + mat[6][5] = 0.0; + mat[6][6] = 1.50*(zz-1.0)*zz; + } + } +} +!------------------------------------------------------------------------------! + SUBROUTINE set_dmat(dmat,r,num) + IMPLICIT NONE + REAL (dbl), INTENT (IN) :: r(3) + INTEGER, INTENT (IN) :: num + REAL (dbl), INTENT (OUT) :: dmat(0:,0:,1:) + INTEGER :: l1, l2 + REAL (dbl), PARAMETER :: s3 = 0.5773502691896260 ! sqrt(1.0/3.0) + REAL (dbl) :: x, y, z, xx, yy, zz, xy, xz, yz + REAL (dbl) :: dc(0:6,0:6) + + l1 = size(dmat(:,0,1)) + l2 = size(dmat(0,:,1)) + x = r(1) + y = r(2) + z = r(3) + IF (l1 == 5 && l2 == 3 && num == 1) THEN +! d1(d,p) + dmat(0:4,0:2,1:3) = zero + dmat(1,0,1) = one + dmat(0,1,1) = -s3 + dmat(3,1,1) = one + dmat(4,2,1) = one + dmat(2,0,2) = one + dmat(0,2,2) = -s3 + dmat(3,2,2) = -one + dmat(4,1,2) = one + dmat(0,0,3) = 2.0*s3 + dmat(1,1,3) = one + dmat(2,2,3) = one + ELSE IF (l1 == 3 && l2 == 5 && num == 1) THEN +! d1(p,d) + dmat(0:2,0:4,1:3) = zero + dmat(0,1,1) = one + dmat(1,0,1) = -s3 + dmat(1,3,1) = one + dmat(2,4,1) = one + dmat(0,2,2) = one + dmat(2,0,2) = -s3 + dmat(2,3,2) = -one + dmat(1,4,2) = one + dmat(0,0,3) = 2.0*s3 + dmat(1,1,3) = one + dmat(2,2,3) = one + ELSEIF (l1 == 5 && l2 == 5 && num == 2) THEN +! d2(d,d) + dmat(0:4,0:4,1:3) = zero + dmat(0,1,1) = s3*z + dmat(0,3,1) = -2.0*s3*x + dmat(0,4,1) = -2.0*s3*y + dmat(1,0,1) = dmat(0,1,1) + dmat(3,0,1) = dmat(0,3,1) + dmat(4,0,1) = dmat(0,4,1) + dmat(1,2,1) = y + dmat(1,3,1) = z + dmat(2,1,1) = y + dmat(3,1,1) = z + dmat(2,2,1) = -2.0*x + dmat(2,4,1) = z + dmat(4,2,1) = z + dmat(0,2,2) = s3*z + dmat(0,3,2) = 2.0*s3*y + dmat(0,4,2) = -2.0*s3*x + dmat(2,0,2) = dmat(0,2,2) + dmat(3,0,2) = dmat(0,3,2) + dmat(4,0,2) = dmat(0,4,2) + dmat(1,1,2) = -2.0*y + dmat(1,2,2) = x + dmat(1,4,2) = z + dmat(2,1,2) = x + dmat(4,1,2) = z + dmat(2,3,2) = -z + dmat(3,2,2) = -z + dmat(0,0,3) = 2.0*z + dmat(0,1,3) = s3*x + dmat(0,2,3) = s3*y + dmat(1,0,3) = dmat(0,1,3) + dmat(2,0,3) = dmat(0,2,3) + dmat(1,3,3) = x + dmat(1,4,3) = y + dmat(3,1,3) = x + dmat(4,1,3) = y + dmat(2,3,3) = -y + dmat(2,4,3) = x + dmat(3,2,3) = -y + dmat(4,2,3) = x + dmat(3,3,3) = -2.0*z + dmat(4,4,3) = -2.0*z + ELSE IF (l1 == 7 && l2 == 3 && num == 2) THEN +! d2(f,p) + dmat(0,0,1) = zero + dmat(0,1,1) = -sqrt(6.0)*z + dmat(0,2,1) = zero + dmat(1,0,1) = 4.0*z + dmat(1,1,1) = -2.0*x + dmat(1,2,1) = -y + dmat(2,0,1) = zero + dmat(2,1,1) = -y + dmat(2,2,1) = zero + dmat(3,0,1) = 2.0*sqrt(2.50)*x + dmat(3,1,1) = sqrt(10.0)*z + dmat(3,2,1) = zero + dmat(4,0,1) = sqrt(10.0)*y + dmat(4,1,1) = zero + dmat(4,2,1) = sqrt(10.0)*z + dmat(5,0,1) = zero + dmat(5,1,1) = sqrt(15.0)*x + dmat(5,2,1) = -sqrt(15.0)*y + dmat(6,0,1) = zero + dmat(6,1,1) = sqrt(15.0)*y + dmat(6,2,1) = sqrt(15.0)*x + + dmat(0,0,2) = zero + dmat(0,1,2) = zero + dmat(0,2,2) = -sqrt(6.0)*z + dmat(1,0,2) = zero + dmat(1,1,2) = zero + dmat(1,2,2) = -x + dmat(2,0,2) = 4.0*z + dmat(2,1,2) = -x + dmat(2,2,2) = -2.0*y + dmat(3,0,2) = -2.0*sqrt(2.50)*y + dmat(3,1,2) = zero + dmat(3,2,2) = -sqrt(10.0)*z + dmat(4,0,2) = sqrt(10.0)*x + dmat(4,1,2) = sqrt(10.0)*z + dmat(4,2,2) = zero + dmat(5,0,2) = zero + dmat(5,1,2) = -sqrt(15.0)*y + dmat(5,2,2) = -sqrt(15.0)*x + dmat(6,0,2) = zero + dmat(6,1,2) = sqrt(15.0)*x + dmat(6,2,2) = -sqrt(15.0)*y + + dmat(0,0,3) = 6.0*sqrt(1.50)*z + dmat(0,1,3) = -sqrt(6.0)*x + dmat(0,2,3) = -sqrt(6.0)*y + dmat(1,0,3) = 4.0*x + dmat(1,1,3) = 5.0*z + dmat(1,2,3) = zero + dmat(2,0,3) = 4.0*y + dmat(2,1,3) = zero + dmat(2,2,3) = 5.0*z + dmat(3,0,3) = zero + dmat(3,1,3) = sqrt(10.0)*x + dmat(3,2,3) = -sqrt(10.0)*y + dmat(4,0,3) = zero + dmat(4,1,3) = sqrt(10.0)*y + dmat(4,2,3) = sqrt(10.0)*x + dmat(5,0,3) = zero + dmat(5,1,3) = zero + dmat(5,2,3) = zero + dmat(6,0,3) = zero + dmat(6,1,3) = zero + dmat(6,2,3) = zero + ELSE IF (l1 == 3 && l2 == 7 && num == 2) THEN +! d2(p,f) + dmat(0,0,1) = zero + dmat(1,0,1) = -sqrt(6.0)*z + dmat(2,0,1) = zero + dmat(0,1,1) = 4.0*z + dmat(1,1,1) = -2.0*x + dmat(2,1,1) = -y + dmat(0,2,1) = zero + dmat(1,2,1) = -y + dmat(2,2,1) = zero + dmat(0,3,1) = 2.0*sqrt(2.50)*x + dmat(1,3,1) = sqrt(10.0)*z + dmat(2,3,1) = zero + dmat(0,4,1) = sqrt(10.0)*y + dmat(1,4,1) = zero + dmat(2,4,1) = sqrt(10.0)*z + dmat(0,5,1) = zero + dmat(1,5,1) = sqrt(15.0)*x + dmat(2,5,1) = -sqrt(15.0)*y + dmat(0,6,1) = zero + dmat(1,6,1) = sqrt(15.0)*y + dmat(2,6,1) = sqrt(15.0)*x + + dmat(0,0,2) = zero + dmat(1,0,2) = zero + dmat(2,0,2) = -sqrt(6.0)*z + dmat(0,1,2) = zero + dmat(1,1,2) = zero + dmat(2,1,2) = -x + dmat(0,2,2) = 4.0*z + dmat(1,2,2) = -x + dmat(2,2,2) = -2.0*y + dmat(0,3,2) = -2.0*sqrt(2.50)*y + dmat(1,3,2) = zero + dmat(2,3,2) = -sqrt(10.0)*z + dmat(0,4,2) = sqrt(10.0)*x + dmat(1,4,2) = sqrt(10.0)*z + dmat(2,4,2) = zero + dmat(0,5,2) = zero + dmat(1,5,2) = -sqrt(15.0)*y + dmat(2,5,2) = -sqrt(15.0)*x + dmat(0,6,2) = zero + dmat(1,6,2) = sqrt(15.0)*x + dmat(2,6,2) = -sqrt(15.0)*y + + dmat(0,0,3) = 6.0*sqrt(1.50)*z + dmat(1,0,3) = -sqrt(6.0)*x + dmat(2,0,3) = -sqrt(6.0)*y + dmat(0,1,3) = 4.0*x + dmat(1,1,3) = 5.0*z + dmat(2,1,3) = zero + dmat(0,2,3) = 4.0*y + dmat(1,2,3) = zero + dmat(2,2,3) = 5.0*z + dmat(0,3,3) = zero + dmat(1,3,3) = sqrt(10.0)*x + dmat(2,3,3) = -sqrt(10.0)*y + dmat(0,4,3) = zero + dmat(1,4,3) = sqrt(10.0)*y + dmat(2,4,3) = sqrt(10.0)*x + dmat(0,5,3) = zero + dmat(1,5,3) = zero + dmat(2,5,3) = zero + dmat(0,6,3) = zero + dmat(1,6,3) = zero + dmat(2,6,3) = zero + ELSE IF (l1 == 7 && l2 == 5 && num == 1) THEN +! d1(f,d) + dmat(0:6,0:4,1:3) = zero + dmat(0,1,1) = -sqrt(1.50) + dmat(1,0,1) = sqrt(3.0) + dmat(1,3,1) = -0.50 + dmat(2,4,1) = -0.50 + dmat(3,1,1) = sqrt(2.50) + dmat(4,2,1) = sqrt(2.50) + dmat(5,3,1) = sqrt(3.750) + dmat(6,4,1) = sqrt(3.750) + + dmat(0,2,2) = -sqrt(1.50) + dmat(1,4,2) = -0.50 + dmat(2,0,2) = sqrt(3.0) + dmat(2,3,2) = 0.50 + dmat(3,2,2) = -sqrt(2.50) + dmat(4,1,2) = sqrt(2.50) + dmat(5,4,2) = -sqrt(3.750) + dmat(6,3,2) = sqrt(3.750) + + dmat(0,0,3) = sqrt(4.50) + dmat(1,1,3) = 2.0 + dmat(2,2,3) = 2.0 + dmat(3,3,3) = sqrt(2.50) + dmat(4,4,3) = sqrt(2.50) + ELSE IF (l1 == 5 && l2 == 7 && num == 1) THEN +! d1(d,f) + dmat(0:4,0:6,1:3) = zero + dmat(1,0,1) = -sqrt(1.50) + dmat(0,1,1) = sqrt(3.0) + dmat(3,1,1) = -0.50 + dmat(4,2,1) = -0.50 + dmat(1,3,1) = sqrt(2.50) + dmat(2,4,1) = sqrt(2.50) + dmat(3,5,1) = sqrt(3.750) + dmat(4,6,1) = sqrt(3.750) + + dmat(2,0,2) = -sqrt(1.50) + dmat(4,1,2) = -0.50 + dmat(0,2,2) = sqrt(3.0) + dmat(3,2,2) = 0.50 + dmat(2,3,2) = -sqrt(2.50) + dmat(1,4,2) = sqrt(2.50) + dmat(4,5,2) = -sqrt(3.750) + dmat(3,6,2) = sqrt(3.750) + + dmat(0,0,3) = sqrt(4.50) + dmat(1,1,3) = 2.0 + dmat(2,2,3) = 2.0 + dmat(3,3,3) = sqrt(2.50) + dmat(4,4,3) = sqrt(2.50) + ELSE IF (l1 == 7 && l2 == 5 && num == 3) THEN + xx = x*x + yy = y*y + zz = z*z + xy = x*y + xz = x*z + yz = y*z +! d3(f,d) + dmat(0,0,1) = zero + dmat(0,1,1) = sqrt(1.50)*zz + dmat(0,2,1) = zero + dmat(0,3,1) = -2.0*sqrt(6.0)*xz + dmat(0,4,1) = -2.0*sqrt(6.0)*yz + dmat(1,0,1) = sqrt(0.750)*(3.0*zz-1.0) + dmat(1,1,1) = xz + dmat(1,2,1) = 3.0*yz + dmat(1,3,1) = 2.50*zz-3.0*xx+yy + dmat(1,4,1) = -4.0*xy + dmat(2,0,1) = zero + dmat(2,1,1) = 3.0*yz + dmat(2,2,1) = -5.0*xz + dmat(2,3,1) = -2.0*xy + dmat(2,4,1) = 0.50*(5.0*zz-4.0*yy) + dmat(3,0,1) = zero + dmat(3,1,1) = sqrt(2.50)*(zz-2.0*yy) + dmat(3,2,1) = sqrt(2.50)*4.0*xy + dmat(3,3,1) = zero + dmat(3,4,1) = zero + dmat(4,0,1) = zero + dmat(4,1,1) = zero + dmat(4,2,1) = sqrt(2.50)*(1.0-6.0*xx) + dmat(4,3,1) = zero + dmat(4,4,1) = zero + dmat(5,0,1) = 3.0*sqrt(1.250)*(yy-xx) + dmat(5,1,1) = sqrt(15.0)*xz + dmat(5,2,1) = -sqrt(15.0)*yz + dmat(5,3,1) = -sqrt(3.750)*zz + dmat(5,4,1) = zero + dmat(6,0,1) = -sqrt(5.0)*3.0*xy + dmat(6,1,1) = sqrt(15.0)*yz + dmat(6,2,1) = sqrt(15.0)*xz + dmat(6,3,1) = zero + dmat(6,4,1) = -sqrt(3.750)*zz + + dmat(0,0,2) = zero + dmat(0,1,2) = zero + dmat(0,2,2) = sqrt(1.50)*zz + dmat(0,3,2) = 2.0*sqrt(6.0)*yz + dmat(0,4,2) = -2.0*sqrt(6.0)*xz + dmat(1,0,2) = zero + dmat(1,1,2) = -5.0*yz + dmat(1,2,2) = 3.0*xz + dmat(1,3,2) = 2.0*xy + dmat(1,4,2) = 0.50*(5.0*zz-4.0*xx) + dmat(2,0,2) = sqrt(0.750)*(3.0*zz-1.0) + dmat(2,1,2) = 3.0*xz + dmat(2,2,2) = yz + dmat(2,3,2) = -(2.50*zz+xx-3.0*yy) + dmat(2,4,2) = -4.0*xy + dmat(3,0,2) = zero + dmat(3,1,2) = -sqrt(2.50)*4.0*xy + dmat(3,2,2) = -sqrt(2.50)*(zz-2.0*xx) + dmat(3,3,2) = zero + dmat(3,4,2) = zero + dmat(4,0,2) = zero + dmat(4,1,2) = sqrt(2.50)*(1.0-6.0*yy) + dmat(4,2,2) = zero + dmat(4,3,2) = zero + dmat(4,4,2) = zero + dmat(5,0,2) = 3.0*sqrt(5.0)*xy + dmat(5,1,2) = -sqrt(15.0)*yz + dmat(5,2,2) = -sqrt(15.0)*xz + dmat(5,3,2) = zero + dmat(5,4,2) = sqrt(3.750)*zz + dmat(6,0,2) = -sqrt(1.250)*(3.0*xx-3.0*yy) + dmat(6,1,2) = sqrt(15.0)*xz + dmat(6,2,2) = -sqrt(15.0)*yz + dmat(6,3,2) = -sqrt(3.750)*zz + dmat(6,4,2) = zero + + dmat(0,0,3) = sqrt(0.50)*(12.0*zz-3.0) + dmat(0,1,3) = 2.0*sqrt(1.50)*xz + dmat(0,2,3) = 2.0*sqrt(1.50)*yz + dmat(0,3,3) = -sqrt(6.0)*(xx-yy) + dmat(0,4,3) = -2.0*sqrt(6.0)*xy + dmat(1,0,3) = 3.0*sqrt(3.0)*xz + dmat(1,1,3) = 0.50*(xx-5.0*yy) + dmat(1,2,3) = 3.0*xy + dmat(1,3,3) = 5.0*xz + dmat(1,4,3) = 5.0*yz + dmat(2,0,3) = 3.0*sqrt(3.0)*yz + dmat(2,1,3) = 3.0*xy + dmat(2,2,3) = 0.50*(yy-5.0*xx) + dmat(2,3,3) = -5.0*yz + dmat(2,4,3) = 5.0*xz + dmat(3,0,3) = zero + dmat(3,1,3) = 2.0*sqrt(2.50)*xz + dmat(3,2,3) = -2.0*sqrt(2.50)*yz + dmat(3,3,3) = sqrt(2.50)*(1.0-6.0*zz) + dmat(3,4,3) = zero + dmat(4,0,3) = zero + dmat(4,1,3) = zero + dmat(4,2,3) = zero + dmat(4,3,3) = zero + dmat(4,4,3) = sqrt(2.50)*(1.0-6.0*zz) + dmat(5,0,3) = zero + dmat(5,1,3) = sqrt(3.750)*(xx-yy) + dmat(5,2,3) = -sqrt(15.0)*xy + dmat(5,3,3) = -sqrt(15.0)*xz + dmat(5,4,3) = sqrt(15.0)*yz + dmat(6,0,3) = zero + dmat(6,1,3) = sqrt(15.0)*xy + dmat(6,2,3) = sqrt(3.750)*(xx-yy) + dmat(6,3,3) = -sqrt(15.0)*yz + dmat(6,4,3) = -sqrt(15.0)*xz + + ELSE IF (l1 == 5 && l2 == 7 && num == 3) THEN + xx = x*x + yy = y*y + zz = z*z + xy = x*y + xz = x*z + yz = y*z +! d3(d,f) + dmat(0,0,1) = zero + dmat(1,0,1) = sqrt(1.50)*zz + dmat(2,0,1) = zero + dmat(3,0,1) = -2.0*sqrt(6.0)*xz + dmat(4,0,1) = -2.0*sqrt(6.0)*yz + dmat(0,1,1) = sqrt(0.750)*(3.0*zz-1.0) + dmat(1,1,1) = xz + dmat(2,1,1) = 3.0*yz + dmat(3,1,1) = 2.50*zz-3.0*xx+yy + dmat(4,1,1) = -4.0*xy + dmat(0,2,1) = zero + dmat(1,2,1) = 3.0*yz + dmat(2,2,1) = -5.0*xz + dmat(3,2,1) = -2.0*xy + dmat(4,2,1) = 0.50*(5.0*zz-4.0*yy) + dmat(0,3,1) = zero + dmat(1,3,1) = sqrt(2.50)*(zz-2.0*yy) + dmat(2,3,1) = sqrt(2.50)*4.0*xy + dmat(3,3,1) = zero + dmat(4,3,1) = zero + dmat(0,4,1) = zero + dmat(1,4,1) = zero + dmat(2,4,1) = sqrt(2.50)*(1.0-6.0*xx) + dmat(3,4,1) = zero + dmat(4,4,1) = zero + dmat(0,5,1) = 3.0*sqrt(1.250)*(yy-xx) + dmat(1,5,1) = sqrt(15.0)*xz + dmat(2,5,1) = -sqrt(15.0)*yz + dmat(3,5,1) = -sqrt(3.750)*zz + dmat(4,5,1) = zero + dmat(0,6,1) = -sqrt(5.0)*3.0*xy + dmat(1,6,1) = sqrt(15.0)*yz + dmat(2,6,1) = sqrt(15.0)*xz + dmat(3,6,1) = zero + dmat(4,6,1) = -sqrt(3.750)*zz + + dmat(0,0,2) = zero + dmat(1,0,2) = zero + dmat(2,0,2) = sqrt(1.50)*zz + dmat(3,0,2) = 2.0*sqrt(6.0)*yz + dmat(4,0,2) = -2.0*sqrt(6.0)*xz + dmat(0,1,2) = zero + dmat(1,1,2) = -5.0*yz + dmat(2,1,2) = 3.0*xz + dmat(3,1,2) = 2.0*xy + dmat(4,1,2) = 0.50*(5.0*zz-4.0*xx) + dmat(0,2,2) = sqrt(0.750)*(3.0*zz-1.0) + dmat(1,2,2) = 3.0*xz + dmat(2,2,2) = yz + dmat(3,2,2) = -(2.50*zz+xx-3.0*yy) + dmat(4,2,2) = -4.0*xy + dmat(0,3,2) = zero + dmat(1,3,2) = -sqrt(2.50)*4.0*xy + dmat(2,3,2) = -sqrt(2.50)*(zz-2.0*xx) + dmat(3,3,2) = zero + dmat(4,3,2) = zero + dmat(0,4,2) = zero + dmat(1,4,2) = sqrt(2.50)*(1.0-6.0*yy) + dmat(2,4,2) = zero + dmat(3,4,2) = zero + dmat(4,4,2) = zero + dmat(0,5,2) = sqrt(5.0)*3.0*xy + dmat(1,5,2) = -sqrt(15.0)*yz + dmat(2,5,2) = -sqrt(15.0)*xz + dmat(3,5,2) = zero + dmat(4,5,2) = sqrt(3.750)*zz + dmat(0,6,2) = -sqrt(1.250)*(3.0*xx-3.0*yy) + dmat(1,6,2) = sqrt(15.0)*xz + dmat(2,6,2) = -sqrt(15.0)*yz + dmat(3,6,2) = -sqrt(3.750)*zz + dmat(4,6,2) = zero + + dmat(0,0,3) = sqrt(0.50)*(12.0*zz-3.0) + dmat(1,0,3) = 2.0*sqrt(1.50)*xz + dmat(2,0,3) = 2.0*sqrt(1.50)*yz + dmat(3,0,3) = -sqrt(6.0)*(xx-yy) + dmat(4,0,3) = -2.0*sqrt(6.0)*xy + dmat(0,1,3) = 3.0*sqrt(3.0)*xz + dmat(1,1,3) = 0.50*(xx-5.0*yy) + dmat(2,1,3) = 3.0*xy + dmat(3,1,3) = 5.0*xz + dmat(4,1,3) = 5.0*yz + dmat(0,2,3) = 3.0*sqrt(3.0)*yz + dmat(1,2,3) = 3.0*xy + dmat(2,2,3) = 0.50*(yy-5.0*xx) + dmat(3,2,3) = -5.0*yz + dmat(4,2,3) = 5.0*xz + dmat(0,3,3) = zero + dmat(1,3,3) = 2.0*sqrt(2.50)*xz + dmat(2,3,3) = -2.0*sqrt(2.50)*yz + dmat(3,3,3) = sqrt(2.50)*(1.0-6.0*zz) + dmat(4,3,3) = zero + dmat(0,4,3) = zero + dmat(1,4,3) = zero + dmat(2,4,3) = zero + dmat(3,4,3) = zero + dmat(4,4,3) = sqrt(2.50)*(1.0-6.0*zz) + dmat(0,5,3) = zero + dmat(1,5,3) = sqrt(3.750)*(xx-yy) + dmat(2,5,3) = -sqrt(15.0)*xy + dmat(3,5,3) = -sqrt(15.0)*xz + dmat(4,5,3) = sqrt(15.0)*yz + dmat(0,6,3) = zero + dmat(1,6,3) = sqrt(15.0)*xy + dmat(2,6,3) = sqrt(3.750)*(xx-yy) + dmat(3,6,3) = -sqrt(15.0)*yz + dmat(4,6,3) = -sqrt(15.0)*xz + ELSE IF (l1 == 7 && l2 == 7 && num == 2) THEN +! d2(f,f) + dmat(0,0,1) = zero + dmat(0,1,1) = sqrt(0.240)*z + dmat(0,2,1) = zero + dmat(0,3,1) = -2.0*sqrt(0.60)*x + dmat(0,4,1) = -sqrt(2.40)*y + dmat(0,5,1) = zero + dmat(0,6,1) = zero + dmat(1,0,1) = dmat(0,1,1) + dmat(2,0,1) = dmat(0,2,1) + dmat(3,0,1) = dmat(0,3,1) + dmat(4,0,1) = dmat(0,4,1) + dmat(5,0,1) = dmat(0,5,1) + dmat(6,0,1) = dmat(0,6,1) + dmat(1,1,1) = zero + dmat(1,2,1) = 1.20*y + dmat(1,3,1) = sqrt(0.90)*z + dmat(1,4,1) = zero + dmat(1,5,1) = -2.0*sqrt(0.150)*x + dmat(1,6,1) = -sqrt(0.60)*y + dmat(2,1,1) = dmat(1,2,1) + dmat(3,1,1) = dmat(1,3,1) + dmat(4,1,1) = dmat(1,4,1) + dmat(5,1,1) = dmat(1,5,1) + dmat(6,1,1) = dmat(1,6,1) + dmat(2,2,1) = -2.40*x + dmat(2,3,1) = zero + dmat(2,4,1) = sqrt(0.90)*z + dmat(2,5,1) = sqrt(0.60)*y + dmat(2,6,1) = -2.0*sqrt(0.150)*x + dmat(3,2,1) = dmat(2,3,1) + dmat(4,2,1) = dmat(2,4,1) + dmat(5,2,1) = dmat(2,5,1) + dmat(6,2,1) = dmat(2,6,1) + dmat(3,3,1) = zero + dmat(3,4,1) = zero + dmat(3,5,1) = sqrt(1.50)*z + dmat(3,6,1) = zero + dmat(4,3,1) = dmat(3,4,1) + dmat(5,3,1) = dmat(3,5,1) + dmat(6,3,1) = dmat(3,6,1) + dmat(4,4,1) = zero + dmat(4,5,1) = zero + dmat(4,6,1) = sqrt(1.50)*z + dmat(5,4,1) = dmat(4,5,1) + dmat(6,4,1) = dmat(4,6,1) + dmat(5,5,1) = zero + dmat(5,6,1) = zero + dmat(6,5,1) = zero + dmat(6,6,1) = zero + + dmat(0,0,2) = zero + dmat(0,1,2) = zero + dmat(0,2,2) = sqrt(0.240)*z + dmat(0,3,2) = 2.0*sqrt(0.60)*y + dmat(0,4,2) = -sqrt(2.40)*x + dmat(0,5,2) = zero + dmat(0,6,2) = zero + dmat(1,0,2) = dmat(0,1,2) + dmat(2,0,2) = dmat(0,2,2) + dmat(3,0,2) = dmat(0,3,2) + dmat(4,0,2) = dmat(0,4,2) + dmat(5,0,2) = dmat(0,5,2) + dmat(6,0,2) = dmat(0,6,2) + dmat(1,1,2) = -2.40*y + dmat(1,2,2) = 1.20*x + dmat(1,3,2) = zero + dmat(1,4,2) = sqrt(0.90)*z + dmat(1,5,2) = 2.0*sqrt(0.150)*y + dmat(1,6,2) = -sqrt(0.60)*x + dmat(2,1,2) = dmat(1,2,2) + dmat(3,1,2) = dmat(1,3,2) + dmat(4,1,2) = dmat(1,4,2) + dmat(5,1,2) = dmat(1,5,2) + dmat(6,1,2) = dmat(1,6,2) + dmat(2,2,2) = zero + dmat(2,3,2) = -sqrt(0.90)*z + dmat(2,4,2) = zero + dmat(2,5,2) = sqrt(0.60)*x + dmat(2,6,2) = 2.0*sqrt(0.150)*y + dmat(3,2,2) = dmat(2,3,2) + dmat(4,2,2) = dmat(2,4,2) + dmat(5,2,2) = dmat(2,5,2) + dmat(6,2,2) = dmat(2,6,2) + dmat(3,3,2) = zero + dmat(3,4,2) = zero + dmat(3,5,2) = zero + dmat(3,6,2) = sqrt(1.50)*z + dmat(4,3,2) = dmat(3,4,2) + dmat(5,3,2) = dmat(3,5,2) + dmat(6,3,2) = dmat(3,6,2) + dmat(4,4,2) = zero + dmat(4,5,2) = -sqrt(1.50)*z + dmat(4,6,2) = zero + dmat(5,4,2) = dmat(4,5,2) + dmat(6,4,2) = dmat(4,6,2) + dmat(5,5,2) = zero + dmat(5,6,2) = zero + dmat(6,5,2) = zero + dmat(6,6,2) = zero + + dmat(0,0,3) = 2.40*z + dmat(0,1,3) = sqrt(0.240)*x + dmat(0,2,3) = sqrt(0.240)*y + dmat(0,3,3) = zero + dmat(0,4,3) = zero + dmat(0,5,3) = zero + dmat(0,6,3) = zero + dmat(1,0,3) = dmat(0,1,3) + dmat(2,0,3) = dmat(0,2,3) + dmat(3,0,3) = dmat(0,3,3) + dmat(4,0,3) = dmat(0,4,3) + dmat(5,0,3) = dmat(0,5,3) + dmat(6,0,3) = dmat(0,6,3) + dmat(1,1,3) = 0.60*z + dmat(1,2,3) = zero + dmat(1,3,3) = sqrt(0.90)*x + dmat(1,4,3) = sqrt(0.90)*y + dmat(1,5,3) = zero + dmat(1,6,3) = zero + dmat(2,1,3) = dmat(1,2,3) + dmat(3,1,3) = dmat(1,3,3) + dmat(4,1,3) = dmat(1,4,3) + dmat(5,1,3) = dmat(1,5,3) + dmat(6,1,3) = dmat(1,6,3) + dmat(2,2,3) = 0.60*z + dmat(2,3,3) = -sqrt(0.90)*y + dmat(2,4,3) = sqrt(0.90)*x + dmat(2,5,3) = zero + dmat(2,6,3) = zero + dmat(3,2,3) = dmat(2,3,3) + dmat(4,2,3) = dmat(2,4,3) + dmat(5,2,3) = dmat(2,5,3) + dmat(6,2,3) = dmat(2,6,3) + dmat(3,3,3) = zero + dmat(3,4,3) = zero + dmat(3,5,3) = sqrt(1.50)*x + dmat(3,6,3) = sqrt(1.50)*y + dmat(4,3,3) = dmat(3,4,3) + dmat(5,3,3) = dmat(3,5,3) + dmat(6,3,3) = dmat(3,6,3) + dmat(4,4,3) = zero + dmat(4,5,3) = -sqrt(1.50)*y + dmat(4,6,3) = sqrt(1.50)*x + dmat(5,4,3) = dmat(4,5,3) + dmat(6,4,3) = dmat(4,6,3) + dmat(5,5,3) = -3.0*z + dmat(5,6,3) = zero + dmat(6,5,3) = zero + dmat(6,6,3) = -3.0*z + ELSE IF (l1 == 7 && l2 == 7 && num == 4) THEN + xx = x*x + yy = y*y + zz = z*z + xy = x*y + xz = x*z + yz = y*z +! d4(f,f) + dmat(0,0,1) = zero + dmat(0,1,1) = sqrt(0.240)*(5.0*zz-2.0)*z + dmat(0,2,1) = zero + dmat(0,3,1) = -sqrt(0.60)*2.0*x*zz + dmat(0,4,1) = -sqrt(2.40)*y*zz + dmat(0,5,1) = sqrt(3.60)*(3.0*yy-3.0*xx)*z + dmat(0,6,1) = -sqrt(3.60)*6.0*x*yz + dmat(1,0,1) = dmat(0,1,1) + dmat(2,0,1) = dmat(0,2,1) + dmat(3,0,1) = dmat(0,3,1) + dmat(4,0,1) = dmat(0,4,1) + dmat(5,0,1) = dmat(0,5,1) + dmat(6,0,1) = dmat(0,6,1) + dmat(1,1,1) = -0.80*x + 3.0*x*zz + dmat(1,2,1) = 0.40*(10.0*zz-1.0)*y + dmat(1,3,1) = sqrt(0.10)*(6.0*zz-8.0*yy-1.0)*z + dmat(1,4,1) = zero + dmat(1,5,1) = sqrt(3.750)*(2.0*x-5.60*x*xx+2.40*x*yy) + dmat(1,6,1) = sqrt(0.60)*(4.0*zz-12.0*xx+1.0)*y + dmat(2,1,1) = dmat(1,2,1) + dmat(3,1,1) = dmat(1,3,1) + dmat(4,1,1) = dmat(1,4,1) + dmat(5,1,1) = dmat(1,5,1) + dmat(6,1,1) = dmat(1,6,1) + dmat(2,2,1) = -5.0*x*zz + dmat(2,3,1) = sqrt(0.10)*16.0*x*yz + dmat(2,4,1) = sqrt(0.10)*(2.0*zz-24.0*xx+3.0)*z + dmat(2,5,1) = sqrt(0.60)*(4.0*yy-4.0*zz-1.0)*y + dmat(2,6,1) = sqrt(3.750)*(2.0*x-4.0*x*xx-2.40*x*yy) + dmat(3,2,1) = dmat(2,3,1) + dmat(4,2,1) = dmat(2,4,1) + dmat(5,2,1) = dmat(2,5,1) + dmat(6,2,1) = dmat(2,6,1) + dmat(3,3,1) = -8.0*x*yy + dmat(3,4,1) = 2.0*(3.0*xx-yy)*y + dmat(3,5,1) = -sqrt(1.50)*(2.0*zz-1.0)*z + dmat(3,6,1) = zero + dmat(4,3,1) = dmat(3,4,1) + dmat(5,3,1) = dmat(3,5,1) + dmat(6,3,1) = dmat(3,6,1) + dmat(4,4,1) = 8.0*x*yy + dmat(4,5,1) = zero + dmat(4,6,1) = -sqrt(1.50)*(2.0*zz-1.0)*z + dmat(5,4,1) = dmat(4,5,1) + dmat(6,4,1) = dmat(4,6,1) + dmat(5,5,1) = zero + dmat(5,6,1) = zero + dmat(6,5,1) = zero + dmat(6,6,1) = zero + + dmat(0,0,2) = zero + dmat(0,1,2) = zero + dmat(0,2,2) = sqrt(0.240)*(5.0*zz-2.0)*z + dmat(0,3,2) = sqrt(0.60)*2.0*y*zz + dmat(0,4,2) = -sqrt(2.40)*x*zz + dmat(0,5,2) = sqrt(3.60)*6.0*y*xz + dmat(0,6,2) = sqrt(3.60)*(3.0*yy-3.0*xx)*z + dmat(1,0,2) = dmat(0,1,2) + dmat(2,0,2) = dmat(0,2,2) + dmat(3,0,2) = dmat(0,3,2) + dmat(4,0,2) = dmat(0,4,2) + dmat(5,0,2) = dmat(0,5,2) + dmat(6,0,2) = dmat(0,6,2) + dmat(1,1,2) = -5.0*y*zz + dmat(1,2,2) = 0.40*(10.0*zz-1.0)*x + dmat(1,3,2) = sqrt(0.10)*(-16.0*y)*xz + dmat(1,4,2) = sqrt(0.10)*(2.0*zz-24.0*yy+3.0)*z + dmat(1,5,2) = sqrt(3.750)*(-2.0*y+2.40*xx*y+4.0*y*yy) + dmat(1,6,2) = sqrt(0.60)*(4.0*zz-4.0*xx+1.0)*x + dmat(2,1,2) = dmat(1,2,2) + dmat(3,1,2) = dmat(1,3,2) + dmat(4,1,2) = dmat(1,4,2) + dmat(5,1,2) = dmat(1,5,2) + dmat(6,1,2) = dmat(1,6,2) + dmat(2,2,2) = -0.80*y + 3.0*y*zz + dmat(2,3,2) = -sqrt(0.10)*(6.0*zz-8.0*xx-1.0)*z + dmat(2,4,2) = zero + dmat(2,5,2) = sqrt(0.60)*(12.0*yy-4.0*zz-1.0)*x + dmat(2,6,2) = sqrt(3.750)*(-2.0*y-2.40*xx*y+5.60*y*yy) + dmat(3,2,2) = dmat(2,3,2) + dmat(4,2,2) = dmat(2,4,2) + dmat(5,2,2) = dmat(2,5,2) + dmat(6,2,2) = dmat(2,6,2) + dmat(3,3,2) = -8.0*xx*y + dmat(3,4,2) = 2.0*(xx-3.0*yy)*x + dmat(3,5,2) = zero + dmat(3,6,2) = -sqrt(1.50)*(2.0*zz-1.0)*z + dmat(4,3,2) = dmat(3,4,2) + dmat(5,3,2) = dmat(3,5,2) + dmat(6,3,2) = dmat(3,6,2) + dmat(4,4,2) = 8.0*xx*y + dmat(4,5,2) = sqrt(1.50)*(2.0*zz-1.0)*z + dmat(4,6,2) = zero + dmat(5,4,2) = dmat(4,5,2) + dmat(6,4,2) = dmat(4,6,2) + dmat(5,5,2) = zero + dmat(5,6,2) = zero + dmat(6,5,2) = zero + dmat(6,6,2) = zero + + dmat(0,0,3) = 0.60*(20.0*zz-8.0)*z + dmat(0,1,3) = sqrt(0.240)*(15.0*zz-2.0)*x + dmat(0,2,3) = sqrt(0.240)*(15.0*zz-2.0)*y + dmat(0,3,3) = -sqrt(0.60)*(xx-yy)*2.0*z + dmat(0,4,3) = -sqrt(2.40)*xy*2.0*z + dmat(0,5,3) = sqrt(3.60)*(3.0*yy-xx)*x + dmat(0,6,3) = sqrt(3.60)*(yy-3.0*xx)*y + dmat(1,0,3) = dmat(0,1,3) + dmat(2,0,3) = dmat(0,2,3) + dmat(3,0,3) = dmat(0,3,3) + dmat(4,0,3) = dmat(0,4,3) + dmat(5,0,3) = dmat(0,5,3) + dmat(6,0,3) = dmat(0,6,3) + dmat(1,1,3) = -(5.0*yy-3.0*xx)*z + dmat(1,2,3) = 8.0*z*xy + dmat(1,3,3) = sqrt(0.10)*(18.0*zz-8.0*yy-1.0)*x + dmat(1,4,3) = sqrt(0.10)*(6.0*zz-8.0*yy+3.0)*y + dmat(1,5,3) = zero + dmat(1,6,3) = sqrt(0.60)*8.0*z*xy + dmat(2,1,3) = dmat(1,2,3) + dmat(3,1,3) = dmat(1,3,3) + dmat(4,1,3) = dmat(1,4,3) + dmat(5,1,3) = dmat(1,5,3) + dmat(6,1,3) = dmat(1,6,3) + dmat(2,2,3) = -(5.0*xx-3.0*yy)*z + dmat(2,3,3) = -sqrt(0.10)*(18.0*zz-8.0*xx-1.0)*y + dmat(2,4,3) = sqrt(0.10)*(6.0*zz-8.0*xx+3.0)*x + dmat(2,5,3) = sqrt(0.60)*(-8.0*z)*xy + dmat(2,6,3) = zero + dmat(3,2,3) = dmat(2,3,3) + dmat(4,2,3) = dmat(2,4,3) + dmat(5,2,3) = dmat(2,5,3) + dmat(6,2,3) = dmat(2,6,3) + dmat(3,3,3) = (4.0-12.0*zz)*z + dmat(3,4,3) = zero + dmat(3,5,3) = -sqrt(1.50)*(6.0*zz-1.0)*x + dmat(3,6,3) = -sqrt(1.50)*(6.0*zz-1.0)*y + dmat(4,3,3) = dmat(3,4,3) + dmat(5,3,3) = dmat(3,5,3) + dmat(6,3,3) = dmat(3,6,3) + dmat(4,4,3) = -2.0*(2.0*zz-1.0)*4.0*z + dmat(4,5,3) = sqrt(1.50)*(6.0*zz-1.0)*y + dmat(4,6,3) = -sqrt(1.50)*(6.0*zz-1.0)*x + dmat(5,4,3) = dmat(4,5,3) + dmat(6,4,3) = dmat(4,6,3) + dmat(5,5,3) = 1.50*(4.0*zz-2.0)*z + dmat(5,6,3) = zero + dmat(6,5,3) = zero + dmat(6,6,3) = 1.50*(4.0*zz-2.0)*z + END IF + dc(0:l1-1,0:l2-1) = r(1)*dmat(:,:,1) + r(2)*dmat(:,:,2) + & + r(3)*dmat(:,:,3) + dmat(:,:,1) = dmat(:,:,1) - r(1)*dc(0:l1-1,0:l2-1) + dmat(:,:,2) = dmat(:,:,2) - r(2)*dc(0:l1-1,0:l2-1) + dmat(:,:,3) = dmat(:,:,3) - r(3)*dc(0:l1-1,0:l2-1) + END SUBROUTINE set_dmat +!------------------------------------------------------------------------------! + END MODULE slater_koster_matr +!------------------------------------------------------------------------------! diff --git a/src/slater_koster_util.F b/src/slater_koster_util.F deleted file mode 100644 index 8891979..0000000 --- a/src/slater_koster_util.F +++ /dev/null @@ -1,175 +0,0 @@ -!------------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright (C) 2000 CP2K developers group ! -!------------------------------------------------------------------------------! -! Some utility functions for the Slater-Koster Modules -! - MODULE slater_koster_util -!------------------------------------------------------------------------------! - USE kinds, ONLY : dbl -! - IMPLICIT NONE -! - PRIVATE - PUBLIC :: sph, dsph, out_prod, out_dprod, dpro -! -!------------------------------------------------------------------------------! -! - CONTAINS -! -!------------------------------------------------------------------------------! - SUBROUTINE sph(l,r,gsl) -! Spherical Harmonics up to f functions -! There is a factor sqrt(4*PI)/sqrt(2*l+1) omitted - IMPLICIT NONE - INTEGER, INTENT (IN) :: l - REAL (dbl), INTENT (IN) :: r(3) - REAL (dbl), INTENT (OUT) :: gsl(0:6) - - REAL (dbl), PARAMETER :: sq3 = 1.732050807568877_dbl ! sqrt(3._dbl) - REAL (dbl), PARAMETER :: sq15 = 3.872983346207417_dbl ! sqrt(15._dbl) - REAL (dbl), PARAMETER :: sq38 = 0.612372435695794_dbl ! sqrt(0.375_dbl) - REAL (dbl), PARAMETER :: sq58 = 0.790569415042095_dbl ! sqrt(0.625_dbl) - - SELECT CASE (l) - CASE (0) - gsl(0) = 1._dbl - CASE (1) - gsl(0) = r(3) - gsl(1) = r(1) - gsl(2) = r(2) - CASE (2) - gsl(0) = 0.5_dbl*(3._dbl*r(3)*r(3)-1._dbl) - gsl(1) = sq3*r(1)*r(3) - gsl(2) = sq3*r(2)*r(3) - gsl(3) = 0.5_dbl*sq3*(r(1)*r(1)-r(2)*r(2)) - gsl(4) = sq3*r(1)*r(2) - CASE (3) - gsl(0) = 0.5_dbl*(5._dbl*r(3)*r(3)-3._dbl)*r(3) - gsl(1) = sq38*(5._dbl*r(3)*r(3)-1._dbl)*r(1) - gsl(2) = sq38*(5._dbl*r(3)*r(3)-1._dbl)*r(2) - gsl(3) = 0.5_dbl*sq15*(r(1)*r(1)-r(2)*r(2))*r(3) - gsl(4) = sq15*r(1)*r(2)*r(3) - gsl(5) = sq58*(r(1)*r(1)-3._dbl*r(2)*r(2))*r(1) - gsl(6) = sq58*(3._dbl*r(1)*r(1)-r(2)*r(2))*r(2) - END SELECT - - END SUBROUTINE sph -!------------------------------------------------------------------------------! - SUBROUTINE dsph(l,r,dgsl) -! Derivatives of Spherical Harmonics up to f functions -! There is a factor sqrt(4*PI)/sqrt(2*l+1)/R omitted - IMPLICIT NONE - INTEGER, INTENT (IN) :: l - REAL (dbl), INTENT (IN) :: r(3) - REAL (dbl), INTENT (OUT) :: dgsl(0:6,1:3) - - REAL (dbl), PARAMETER :: sq3 = 1.732050807568877_dbl ! sqrt(3._dbl) - REAL (dbl), PARAMETER :: sq15 = 3.872983346207417_dbl ! sqrt(15._dbl) - REAL (dbl), PARAMETER :: sq38 = 0.612372435695794_dbl ! sqrt(0.375_dbl) - REAL (dbl), PARAMETER :: sq58 = 0.790569415042095_dbl ! sqrt(0.625_dbl) - REAL (dbl) :: dc(0:6) - - SELECT CASE (l) - CASE (0) - dgsl(0,:) = 0._dbl - CASE (1) - dgsl(0:2,1:3) = 0._dbl - dgsl(1,1) = 1._dbl - dgsl(2,2) = 1._dbl - dgsl(0,3) = 1._dbl - CASE (2) - dgsl(0,1) = 0._dbl - dgsl(1,1) = sq3*r(3) - dgsl(2,1) = 0._dbl - dgsl(3,1) = sq3*r(1) - dgsl(4,1) = sq3*r(2) - dgsl(0,2) = 0._dbl - dgsl(1,2) = 0._dbl - dgsl(2,2) = sq3*r(3) - dgsl(3,2) = -sq3*r(2) - dgsl(4,2) = sq3*r(1) - dgsl(0,3) = 3._dbl*r(3) - dgsl(1,3) = sq3*r(1) - dgsl(2,3) = sq3*r(2) - dgsl(3,3) = 0._dbl - dgsl(4,3) = 0._dbl - CASE (3) - dgsl(0,1) = 0._dbl - dgsl(1,1) = sq38*(5._dbl*r(3)*r(3)-1._dbl) - dgsl(2,1) = 0._dbl - dgsl(3,1) = sq15*r(1)*r(3) - dgsl(4,1) = sq15*r(2)*r(3) - dgsl(5,1) = 3._dbl*sq58*(r(1)*r(1)-r(2)*r(2)) - dgsl(6,1) = 6._dbl*sq58*r(1)*r(2) - dgsl(0,2) = 0._dbl - dgsl(1,2) = 0._dbl - dgsl(2,2) = sq38*(5._dbl*r(3)*r(3)-1._dbl) - dgsl(3,2) = -sq15*r(2)*r(3) - dgsl(4,2) = sq15*r(1)*r(3) - dgsl(5,2) = -6._dbl*sq58*r(2)*r(1) - dgsl(6,2) = 3._dbl*sq58*(r(1)*r(1)-r(2)*r(2)) - dgsl(0,3) = 7.5_dbl*r(3)*r(3)-1.5_dbl - dgsl(1,3) = 10._dbl*sq38*r(3)*r(1) - dgsl(2,3) = 10._dbl*sq38*r(3)*r(2) - dgsl(3,3) = 0.5_dbl*sq15*(r(1)*r(1)-r(2)*r(2)) - dgsl(4,3) = sq15*r(1)*r(2) - dgsl(5,3) = 0._dbl - dgsl(6,3) = 0._dbl - END SELECT - dc(0:2*l) = r(1)*dgsl(0:2*l,1) + r(2)*dgsl(0:2*l,2) + r(3)*dgsl(0:2*l,3) - dgsl(0:2*l,1) = dgsl(0:2*l,1) - r(1)*dc(0:2*l) - dgsl(0:2*l,2) = dgsl(0:2*l,2) - r(2)*dc(0:2*l) - dgsl(0:2*l,3) = dgsl(0:2*l,3) - r(3)*dc(0:2*l) - - END SUBROUTINE dsph -!------------------------------------------------------------------------------! - SUBROUTINE out_prod(mat,v1,v2) - IMPLICIT NONE - REAL (dbl), INTENT (IN) :: v1(:), v2(:) - REAL (dbl), INTENT (OUT) :: mat(:,:) - INTEGER :: n1, n2, i, j - - n1 = size(v1) - n2 = size(v2) - DO j = 1, n2 - DO i = 1, n1 - mat(i,j) = v1(i)*v2(j) - END DO - END DO - END SUBROUTINE out_prod -!------------------------------------------------------------------------------! - SUBROUTINE out_dprod(dmat,v1,v2,dv1,dv2) - IMPLICIT NONE - REAL (dbl), INTENT (IN) :: v1(:), v2(:) - REAL (dbl), INTENT (IN) :: dv1(:), dv2(:) - REAL (dbl), INTENT (OUT) :: dmat(:,:) - INTEGER :: n1, n2, i, j - - n1 = size(v1) - n2 = size(v2) - DO j = 1, n2 - DO i = 1, n1 - dmat(i,j) = dv1(i)*v2(j) + v1(i)*dv2(j) - END DO - END DO - END SUBROUTINE out_dprod -!------------------------------------------------------------------------------! - FUNCTION dpro(g1,g2,l1,l2) - IMPLICIT NONE - REAL (dbl), INTENT (IN) :: g1(0:6,0:6), g2(0:6,0:6) - INTEGER, INTENT (IN) :: l1, l2 - REAL (dbl) :: dpro - INTEGER :: i, j - - dpro = 0._dbl - DO i = 0, 2*l2 - DO j = 0, 2*l1 - dpro = dpro + g1(j,i)*g2(j,i) - END DO - END DO - - END FUNCTION dpro -!------------------------------------------------------------------------------! - END MODULE slater_koster_util -!------------------------------------------------------------------------------! diff --git a/src/slater_koster_util.cpp b/src/slater_koster_util.cpp new file mode 100644 index 0000000..9da14f9 --- /dev/null +++ b/src/slater_koster_util.cpp @@ -0,0 +1,161 @@ +#include + +/*----------------------------------------------------------------------------*/ +/* CP2K: A general program to perform molecular dynamics simulations */ +/* Copyright (C) 2000 CP2K developers group */ +/*----------------------------------------------------------------------------*/ +// Some utility functions for the Slater-Koster Modules +// + +namespace slater_koster_util { + + // Spherical Harmonics up to f functions + // There is a factor sqrt(4*PI)/sqrt(2*l+1) omitted + void sph(int l, double r[3], double gsl[7]) { + const double sq3 = 1.732050807568877; // sqrt(3) + const double sq15 = 3.872983346207417; // sqrt(15) + const double sq38 = 0.612372435695794; // sqrt(0.375) + const double sq58 = 0.790569415042095; // sqrt(0.625) + + switch (l) { + case 0: + gsl[0] = 1.0; + break; + case 1: + gsl[0] = r[2]; + gsl[1] = r[0]; + gsl[2] = r[1]; + break; + case 2: + gsl[0] = 0.5 * (3.0 * r[2] * r[2] - 1.0); + gsl[1] = sq3 * r[0] * r[2]; + gsl[2] = sq3 * r[1] * r[2]; + gsl[3] = 0.5 * sq3 * (r[0] * r[0] - r[1] * r[1]); + gsl[4] = sq3 * r[0] * r[1]; + break; + case 3: + gsl[0] = 0.5 * (5.0 * r[2] * r[2] - 3.0) * r[2]; + gsl[1] = sq38 * (5.0 * r[2] * r[2] - 1.0) * r[0]; + gsl[2] = sq38 * (5.0 * r[2] * r[2] - 1.0) * r[1]; + gsl[3] = 0.5 * sq15 * (r[0] * r[0] - r[1] * r[1]) * r[2]; + gsl[4] = sq15 * r[0] * r[1] * r[2]; + gsl[5] = sq58 * (r[0] * r[0] - 3.0 * r[1] * r[1]) * r[0]; + gsl[6] = sq58 * (3.0 * r[0] * r[0] - r[1] * r[1]) * r[1]; + break; + } + } + + // Derivatives of Spherical Harmonics up to f functions + // There is a factor sqrt(4*PI)/sqrt(2*l+1)/R omitted + void dsph(int l, double r[3], double dgsl[7][3]) { + const double sq3 = 1.732050807568877; // sqrt(3) + const double sq15 = 3.872983346207417; // sqrt(15) + const double sq38 = 0.612372435695794; // sqrt(0.375) + const double sq58 = 0.790569415042095; // sqrt(0.625) + double dc[7]; + + switch (l) { + case 0: + dgsl[0][0] = 0.0; + dgsl[0][1] = 0.0; + dgsl[0][2] = 0.0; + break; + case 1: + dgsl[0][0] = 0.0; + dgsl[0][1] = 1.0; + dgsl[0][2] = 0.0; + dgsl[1][0] = 0.0; + dgsl[1][1] = 0.0; + dgsl[1][2] = 0.0; + dgsl[2][0] = 0.0; + dgsl[2][1] = 0.0; + dgsl[2][2] = 1.0; + break; + case 2: + dgsl[0][0] = 0.0; + dgsl[0][1] = sq3 * r[2]; + dgsl[0][2] = 0.0; + dgsl[1][0] = 0.0; + dgsl[1][1] = 0.0; + dgsl[1][2] = -sq3 * r[2]; + dgsl[2][0] = sq3 * r[2]; + dgsl[2][1] = sq3 * r[0]; + dgsl[2][2] = sq3 * r[1]; + dgsl[3][0] = sq3 * r[0]; + dgsl[3][1] = -sq3 * r[2]; + dgsl[3][2] = sq3 * r[1]; + dgsl[4][0] = sq3 * r[1]; + dgsl[4][1] = sq3 * r[2]; + dgsl[4][2] = sq3 * r[0]; + break; + case 3: + dgsl[0][0] = 0.0; + dgsl[0][1] = sq38 * (5.0 * r[2] * r[2] - 1.0); + dgsl[0][2] = 0.0; + dgsl[1][0] = 0.0; + dgsl[1][1] = 0.0; + dgsl[1][2] = -sq38 * (5.0 * r[2] * r[2] - 1.0); + dgsl[2][0] = sq38 * (5.0 * r[2] * r[2] - 1.0); + dgsl[2][1] = sq38 * (5.0 * r[0] * r[2]); + dgsl[2][2] = sq38 * (5.0 * r[1] * r[2]); + dgsl[3][0] = sq15 * r[0] * r[2]; + dgsl[3][1] = -sq15 * r[2] * r[2] + sq15 * r[0] * r[0]; + dgsl[3][2] = sq15 * r[1] * r[2]; + dgsl[4][0] = sq15 * r[1] * r[2]; + dgsl[4][1] = sq15 * r[2] * r[2] - sq15 * r[1] * r[1]; + dgsl[4][2] = -sq15 * r[0] * r[2]; + dgsl[5][0] = 3.0 * sq58 * (r[0] * r[0] - r[1] * r[1]); + dgsl[5][1] = -6.0 * sq58 * r[0] * r[1]; + dgsl[5][2] = 0.0; + dgsl[6][0] = 6.0 * sq58 * r[0] * r[1]; + dgsl[6][1] = 3.0 * sq58 * (r[1] * r[1] - r[0] * r[0]); + dgsl[6][2] = 0.0; + break; + } + + for (int i = 0; i <= 2 * l; i++) { + dc[i] = r[0] * dgsl[i][0] + r[1] * dgsl[i][1] + r[2] * dgsl[i][2]; + dgsl[i][0] -= r[0] * dc[i]; + dgsl[i][1] -= r[1] * dc[i]; + dgsl[i][2] -= r[2] * dc[i]; + } + } + + // Compute the outer product of two vectors + void out_prod(double mat[][3], double v1[], double v2[]) { + int n1 = sizeof(v1) / sizeof(v1[0]); + int n2 = sizeof(v2) / sizeof(v2[0]); + + for (int j = 0; j < n2; j++) { + for (int i = 0; i < n1; i++) { + mat[i][j] = v1[i] * v2[j]; + } + } + } + + // Compute the outer product of two vectors with their derivatives + void out_dprod(double dmat[][3], double v1[], double v2[], double dv1[], double dv2[]) { + int n1 = sizeof(v1) / sizeof(v1[0]); + int n2 = sizeof(v2) / sizeof(v2[0]); + + for (int j = 0; j < n2; j++) { + for (int i = 0; i < n1; i++) { + dmat[i][j] = dv1[i] * v2[j] + v1[i] * dv2[j]; + } + } + } + + // Compute the dot product of two matrices + double dpro(double g1[][7], double g2[][7], int l1, int l2) { + double dpro = 0.0; + + for (int i = 0; i <= 2 * l2; i++) { + for (int j = 0; j <= 2 * l1; j++) { + dpro += g1[j][i] * g2[j][i]; + } + } + + return dpro; + } + +} // namespace slater_koster_util