aparler/cp2k-c
1
0
Fork 0
forked from comp-chem/cp2k-c
Code Pull requests Activity

Translated some more, particularly Slater-Koster

Browse source
This commit is contained in:
Adam Parler 2024-07-19 02:59:42 -07:00
parent 4d602323b0
commit 0d74c00a83
14 changed files with 2064 additions and 2209 deletions
Download patch file Download diff file Expand all files Collapse all files

120
src/amoeba.F
View file

@ -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<ftol) THEN
CALL swap(y(1),y(ilo))
CALL swap(p(1,:),p(ilo,:))
RETURN
END IF
IF (iter>=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<y(ihi)) THEN
y(ihi) = ytry
psum(:) = psum(:) - p(ihi,:) + ptry(:)
p(ihi,:) = ptry(:)
END IF
amotry = ytry
END FUNCTION amotry
END SUBROUTINE amoeba_evaluate
!******************************************************************************
END MODULE amoeba

96
src/amoeba.cpp Normal file
View file

@ -0,0 +1,96 @@
/*---------------------------------------------------------------------------*/
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart */
/*---------------------------------------------------------------------------*/
#include <cmath>
#include <vector>
#include "amoeba.h"
namespace amoeba {
/*---------------------------------------------------------------------------*/
/* NUMERICAL RECIPIES SUBROUTINE AMOEBA ADAPTED FOR USE IN CP2K */
/*---------------------------------------------------------------------------*/
template<typename Func>
void amoeba_evaluate(std::vector<std::vector<double>>& p, std::vector<double>& y,
double ftol, double& rtol, Func func, int itmax) {
int ndim = p[0].size();
std::vector<double> 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<double> 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

10
src/amoeba.h Normal file
View file

@ -0,0 +1,10 @@
#ifndef _AMOEBA_H
#define _AMOEBA_H
#include <vector>
template<typename Func>
void amoeba_evaluate(std::vector<std::vector<double>>&, std::vector<double>&,
double, double&, Func, int);
#endif

246
src/atoms_input.F
View file

@ -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

162
src/atoms_input.cpp Normal file
View file

@ -0,0 +1,162 @@
#include <stdio.h>
#include <stdlib.h>
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]);
}
}

159
src/header.F
View file

@ -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

119
src/header.c Normal file
View file

@ -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");
}
/*****************************************************************************/

21
src/header.h Normal file
View file

@ -0,0 +1,21 @@
#ifndef _HEADER_H
#define _HEADER_H
// Secure block
#ifdef __cplusplus
extern "C" {
#endif
#include <stdio.h>
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

1
src/kinds.cpp
View file

@ -13,6 +13,7 @@
- HITACHI
*/
#include <limits>
#include <iostream>
namespace kinds {

36
src/slater_koster.h Normal file
View file

@ -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

1509
src/slater_koster_matr.F
View file

File diff suppressed because it is too large Load diff

1458
src/slater_koster_matr.c Normal file
View file

File diff suppressed because it is too large Load diff

175
src/slater_koster_util.F
View file

@ -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
!------------------------------------------------------------------------------!

161
src/slater_koster_util.cpp Normal file
View file

@ -0,0 +1,161 @@
#include <cmath>
/*----------------------------------------------------------------------------*/
/* 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