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

Updated more files to C. Mainly ewald and associated files.

Browse source
This commit is contained in:
Adam Parler 2024-08-30 02:57:44 -07:00
parent 0d74c00a83
commit 3e59a77ea8
134 changed files with 3203 additions and 5366 deletions
Download patch file Download diff file Expand all files Collapse all files

4
COPYRIGHT
View file

@ -1,7 +1,7 @@
!-----------------------------------------------------------------------------!
! !
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
! Copyright (C) 2024 CP2K developers group !
! !
! This program is free software; you can redistribute it and/or modify !
! it under the terms of the GNU General Public License as published by !
@ -21,4 +21,4 @@
--
cp2k@prr.mpi-stuttgart.mpg.de
15. X AD 2000
15. X AD 2024

6
INSTALL
View file

@ -4,13 +4,13 @@ How to compile CP2K code
1) Extract the archive file:
# gtar zxvf cp2k.tar.gz
# tar xzvf cp2k.tar.gz
or if you do not have GNU tar
# gzip -cd cp2k.tar.gz | tar xvf -
# unzip -cd cp2k.tar.gz
This extract the contents of the archive into the directory ./CP2K; please be
This extracts the contents of the archive into the directory ./CP2K; please be
careful not to overwrite things in there if it already exists
2) Go to 'CP2K/tools' and do a 'make' (or 'gmake'); this compiles you an

2
Makefile
View file

@ -1,4 +1,4 @@
.SUFFIXES: .cpp .d .o
.SUFFIXES: .c .cpp .d .o
CP2KHOME = $(HOME)/CP2K
ARCH = $(shell $(CP2KHOME)/tools/get_arch_code)

18
src/Makefile
View file

@ -1,4 +1,4 @@
.SUFFIXES: .o .cpp .d
.SUFFIXES: .o .c .cpp .d
CP2KHOME= ${HOME}/CP2K
FORPAR = $(CP2KHOME)/tools/forpar.x -chkint
@ -19,19 +19,25 @@ $(PROG): $(OBJECTS)
$(LDR) $(LFLAGS) -o $(PROG) $(OBJECTS) $(LIBS)
%.o : %.cpp
$(FC) -c $(FFLAGS) $*.cpp
$(CPP) -P $(CDEFS) $*.cpp $*.cpp
$(CXX) -c $(CXXFLAGS) $*.cpp
$(CPP) -P $(CDEFS) $*.cpp
-$(FORPAR) $*.cpp
rm -f $*.cpp
%.o : %.c
$(CC) -c $(CFLAGS) $*.c
$(CC) -P $(CDEFS) $*.c
-$(FORPAR) $*.c
rm -f $*.c
parallel_include.o : parallel_include.cpp
$(FCfixed) -c $(FFLAGS) $*.cpp
$(CXXfixed) -c $(CXXFLAGS) $*.cpp
$(CPP) -P $(CDEFS) $*.cpp $*.cpp
-$(FORPAR) -fix $*.cpp
rm -f $*.cpp
%.d : %.cpp
$(CPP) -P $(CDEFS) $*.cpp $*.cpp
$(CPP) -P $(CDEFS) $*.cpp
$(PERL) $(CP2KHOME)/tools/sfmakedepend -m int -s -f $*.d $*.cpp
rm -f $*.cpp $*.d.old
@ -42,7 +48,7 @@ clean:
rm -f *.o *.mod
realclean:
rm -f $(PROG) *.o *.mod *.d *.int *~ F*.cpp *.lst
rm -f $(PROG) *.o *.mod *.d *.int *~ F*.cpp F*.c *.lst
####################

10
src/amoeba.cpp
View file

@ -1,18 +1,20 @@
/*---------------------------------------------------------------------------*/
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 1999 MPI fuer Festkoerperforschung, Stuttgart */
/* Copyright (C) 2024 MPI fuer Festkoerperforschung, Stuttgart */
/*---------------------------------------------------------------------------*/
#include <cmath>
#include <vector>
#include <algorithm>
#include <functional>
#include "amoeba.h"
namespace amoeba {
/*---------------------------------------------------------------------------*/
/* NUMERICAL RECIPIES SUBROUTINE AMOEBA ADAPTED FOR USE IN CP2K */
/*---------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------|
| 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,

307
src/atoms_input.cpp
View file

@ -1,162 +1,181 @@
#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;
#include "atoms_input.h"
#include "global_types.h"
void read_coord_vel(system_type* atype, char* filen) {
int ierror, iw, source, allgrp, ia, ie, i, ios;
char label[6];
/*----------------------------------------------------------------------------|
| CP2K: A general program to perform molecular dynamics simulations |
| Copyright (C) 2024 CP2K developers group |
|----------------------------------------------------------------------------*/
iw = globenv.scr;
/*>---------------------------------------------------------------------------|
| |
| 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 |
| |
| |
| |
|<---------------------------------------------------------------------------*/
//..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*));
void read_coord_vel(system_type *atype, char *filen) {
int ierror, ilen, iw, source, allgrp, ia, ie, i, ios;
char string[7];
char label[6];
iw = globenv.scr;
//..defaults
//..parse the input section
strncpy(label, "&ATOMS", 10);
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->v[i] = (double*)malloc(atype->n * sizeof(double));
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) {
fprintf(iw, "ATOM| Box parameters [Angstrom]\n");
fprintf(iw, "ATOM| %15.5f %15.5f %15.5f\n", atype->box[0][0], atype->box[0][1],
atype->box[0][2]);
fprintf(iw, "ATOM| %15.5f %15.5f %15.5f\n", atype->box[1][0], atype->box[1][1],
atype->box[1][2]);
fprintf(iw, "ATOM| %15.5f %15.5f %15.5f\n", atype->box[2][0], atype->box[2][1],
atype->box[2][2]);
fprintf(iw, "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);
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");
}
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 read_file(char* filen, system_type* atype) {
int iunit, i, j, ios;
int exists;
void print_c(int iw, double **c, int n) {
int i;
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");
}
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_c(int iw, double** c, int n) {
int 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 %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]);
}
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]);
}
}

33
src/atoms_input.h Normal file
View file

@ -0,0 +1,33 @@
#ifndef _ATOMS_INPUT_H
#define _ATOMS_INPUT_H
#ifdef __cplusplus
namespace CP2K_NS {
#include <string>
#include <vector>
#else
#include <stdbool.h>
#include <string.h>
#endif
typedef struct system_type {
double box[3][3];
char ptype[4];
char rtype[21];
int n;
double** c;
double** v;
} system_type;
void read_coord_vel(system_type, char*, global_environment_type);
void read_file(FILE);
void read_file(FILE, bool);
void print_c(int, std::vector<std::vector<double>>);
void print_cv(int, std::vector<std::vector<double>>, std::vector<std::vector<double>>);
#ifdef __cplusplus
}
#endif
#endif

164
src/coefficient_types.F
View file

@ -1,164 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE coefficient_types
USE kinds, ONLY : dbl
USE pw_grid_types, ONLY : pw_grid_type
USE pw_types, ONLY : pw_type, pw_zero, pw_allocate, pw_deallocate, &
REALDATA1D, COMPLEXDATA1D, REALDATA3D, COMPLEXDATA3D
USE stop_program, ONLY : stop_prg
IMPLICIT NONE
PRIVATE
PUBLIC :: coeff_type, coeff_allocate, coeff_deallocate, coeff_zero
PUBLIC :: PLANEWAVES, GAUSSIAN
PUBLIC :: PW_REALDATA1D, PW_COMPLEXDATA1D, PW_REALDATA3D, PW_COMPLEXDATA3D
!! 'group_id' to grid?
TYPE coeff_type
TYPE ( pw_type ) :: pw
! TYPE ( gaussian_type ) :: gaussian
INTEGER :: in_use
END TYPE coeff_type
!! Flags for the structure member 'in_use'
INTEGER, PARAMETER :: PLANEWAVES = 401, GAUSSIAN = 402
INTERFACE coeff_allocate
MODULE PROCEDURE &
coeff_allocate_from_coeff, coeff_allocate_from_grid
END INTERFACE
!! Flags for 'integral'
PUBLIC :: SQUARE, SQUAREROOT
INTEGER, PARAMETER :: SQUARE = 391, SQUAREROOT = 392
INTEGER, PARAMETER :: PW_REALDATA1D = REALDATA1D
INTEGER, PARAMETER :: PW_COMPLEXDATA1D = COMPLEXDATA1D
INTEGER, PARAMETER :: PW_REALDATA3D = REALDATA3D
INTEGER, PARAMETER :: PW_COMPLEXDATA3D = COMPLEXDATA3D
CONTAINS
!******************************************************************************
SUBROUTINE coeff_allocate_from_coeff ( c_out, c_model, use_basis )
IMPLICIT NONE
! Arguments
TYPE ( coeff_type ), INTENT ( IN ), TARGET :: c_model
TYPE ( coeff_type ), INTENT ( OUT ) :: c_out
INTEGER, INTENT ( IN ), OPTIONAL :: use_basis
! Locals
INTEGER :: use_data
!------------------------------------------------------------------------------
IF ( PRESENT ( use_basis ) ) THEN
IF ( use_basis == PLANEWAVES ) THEN
c_out % in_use = PLANEWAVES
CALL pw_allocate ( c_out % pw, c_model % pw % pw_grid, &
use_data = c_model % pw % in_use )
ELSE IF ( use_basis == PW_REALDATA1D .or. use_basis == PW_COMPLEXDATA1D &
.or. use_basis == PW_REALDATA3D .or. use_basis == PW_COMPLEXDATA3D &
) THEN
c_out % in_use = PLANEWAVES
CALL pw_allocate ( c_out % pw, c_model % pw % pw_grid, &
use_data = use_basis )
ELSE
CALL stop_prg ( "coeff_allocate_from_coeff", &
"no suitable data with use_basis" )
END IF
ELSE
use_data = c_model % pw % in_use
IF ( use_data == PW_REALDATA1D .or. use_data == PW_COMPLEXDATA1D .or. &
use_data == PW_REALDATA3D .or. use_data == PW_COMPLEXDATA3D ) THEN
c_out % in_use = PLANEWAVES
CALL pw_allocate ( c_out % pw, c_model % pw % pw_grid, &
use_data = c_model % pw % in_use )
ELSE
CALL stop_prg ( "coeff_allocate_from_coeff", &
"no suitable data without use_basis" )
END IF
END IF
END SUBROUTINE coeff_allocate_from_coeff
!******************************************************************************
SUBROUTINE coeff_allocate_from_grid ( c_out, pw_grid, use_data )
IMPLICIT NONE
! Arguments
TYPE ( pw_grid_type ), INTENT ( IN ), OPTIONAL :: pw_grid
INTEGER, INTENT ( IN ) :: use_data
TYPE ( coeff_type ), INTENT ( OUT ) :: c_out
! Locals
!------------------------------------------------------------------------------
IF ( PRESENT ( pw_grid ) ) THEN
c_out % in_use = PLANEWAVES
CALL pw_allocate ( c_out % pw, pw_grid, use_data )
ELSE
CALL stop_prg ( "coeff_allocate_from_grid", "no suitable grid" )
END IF
END SUBROUTINE coeff_allocate_from_grid
!******************************************************************************
SUBROUTINE coeff_deallocate ( coeff )
IMPLICIT NONE
! Arguments
TYPE ( coeff_type ), INTENT ( INOUT ) :: coeff
!------------------------------------------------------------------------------
IF ( coeff % in_use == PLANEWAVES ) THEN
CALL pw_deallocate ( coeff % pw )
ELSE
CALL stop_prg ( "coeff_deallocate", "no valid data type" )
END IF
END SUBROUTINE coeff_deallocate
!******************************************************************************
SUBROUTINE coeff_zero ( coeff )
IMPLICIT NONE
! Arguments
TYPE ( coeff_type ), INTENT ( INOUT ) :: coeff
!------------------------------------------------------------------------------
IF ( coeff % in_use == PLANEWAVES ) THEN
CALL pw_zero ( coeff % pw )
ELSE
CALL stop_prg ( "coeff_zero", "no valid data type" )
END IF
END SUBROUTINE coeff_zero
!******************************************************************************
END MODULE coefficient_types

60
src/coefficient_types.cpp Normal file
View file

@ -0,0 +1,60 @@
/*---------------------------------------------------------------------------*/
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 2024 CP2K developers group */
/*---------------------------------------------------------------------------*/
#include "coefficients.h"
using namespace CP2K_NS;
void coeff_allocate_from_coeff(coeff_type c_out, coeff_type c_model, int use_basis) {
int use_data;
if (use_basis > 0) {
if (use_basis == PLANEWAVES) {
c_out.in_use = PLANEWAVES;
pw_allocate(c_out.pw, c_model.pw.pw_grid, c_model.pw.in_use);
} else if (use_basis == PW_REALDATA1D || use_basis == PW_COMPLEXDATA1D ||
use_basis == PW_REALDATA3D || use_basis == PW_COMPLEXDATA3D) {
c_out.in_use = PLANEWAVES;
pw_allocate(c_out.pw, c_model.pw.pw_grid, use_basis);
} else {
stop_prg ("coeff_allocate_from_coeff", "no suitable data with use_basis");
}
} else {
use_data = c_model.pw.in_use;
if ( use_data == PW_REALDATA1D || use_data == PW_COMPLEXDATA1D ||
use_data == PW_REALDATA3D || use_data == PW_COMPLEXDATA3D ) {
c_out.in_use = PLANEWAVES;
pw_allocate(c_out.pw, c_model.pw.pw_grid, c_model.pw.in_use);
} else {
stop_prg("coeff_allocate_from_coeff", "no suitable data without use_basis");
}
}
}
void coeff_allocate_from_grid(coeff_type c_out, pw_grid_type pw_grid, int use_data) {
if (pw_grid) {
c_out.in_use = PLANEWAVES;
pw_allocate(c_out.pw, pw_grid, use_data);
} else {
stop_prg("coeff_allocate_from_grid", "no suitable grid");
}
}
void coeff_deallocate(coeff_type coeff) {
if (coeff.in_use == PLANEWAVES) {
pw_deallocate(coeff.pw);
} else {
stop_prg("coeff_deallocate", "no valid data type");
}
}
void coeff_zero(coeff_type coeff) {
if (coeff.in_use == PLANEWAVES) {
pw_zero(coeff.pw);
} else {
stop_prg("coeff_zero", "no valid data type");
}
}

63
src/coefficients.F
View file

@ -1,63 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE coefficients
USE coefficient_types, ONLY : coeff_type, PLANEWAVES
USE kinds, ONLY : dbl
USE pws, ONLY : pw_integral_a, pw_integral_ab
IMPLICIT NONE
PRIVATE
PUBLIC :: integral
INTERFACE integral
MODULE PROCEDURE integral_a, integral_ab
END INTERFACE
CONTAINS
!******************************************************************************
FUNCTION integral_a ( coeff, flag ) RESULT ( integral_value )
TYPE ( coeff_type ), INTENT ( IN ) :: coeff
INTEGER, INTENT ( IN ), OPTIONAL :: flag
REAL ( dbl ) :: integral_value
!------------------------------------------------------------------------------
IF ( coeff % in_use == PLANEWAVES ) THEN
integral_value = pw_integral_a ( coeff % pw, flag )
ELSE
stop "integral_a, not defined"
END IF
END FUNCTION integral_a
!******************************************************************************
FUNCTION integral_ab ( coeff1, coeff2, flag ) RESULT ( integral_value )
TYPE ( coeff_type ), INTENT ( IN ) :: coeff1, coeff2
INTEGER, INTENT ( IN ), OPTIONAL :: flag
REAL ( dbl ) :: integral_value
!------------------------------------------------------------------------------
IF ( coeff1 % in_use == PLANEWAVES .AND. coeff2 % in_use == PLANEWAVES ) THEN
integral_value = pw_integral_ab ( coeff1 % pw, coeff2 % pw, flag )
ELSE
stop "integral_ab, not defined"
END IF
END FUNCTION integral_ab
!******************************************************************************
END MODULE coefficients

28
src/coefficients.cpp Normal file
View file

@ -0,0 +1,28 @@
/*---------------------------------------------------------------------------*/
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 2024 CP2K developers group */
/*---------------------------------------------------------------------------*/
#include "coefficients.h"
double integral_a(coeff_type *coeff, int flag) {
if (strncpy(coeff->in_use, PLANEWAVES, 256)) {
integral_value = pw_integral_a(coeff->pw, flag);
} else {
printf("integral_a, not defined");
exit(1);
}
return integral_value;
}
double integral_ab(coeff_type *coeff1, coeff_type *coeff2, int flag) {
if (strncpy(coeff1->in_use, PLANEWAVES, 256) && strncpy(coeff2->in_use, PLANEWAVES, 256)) {
integral_value = pw_integral_ab(coeff1->pw, coeff2->pw, flag);
} else {
printf("integral_ab, not defined");
exit(1);
}
return integral_value;
}

36
src/coefficients.h Normal file
View file

@ -0,0 +1,36 @@
#ifndef _COEFF_H
#define _COEFF_H
#include "pw.h"
#include "stop.h"
#ifdef __cplusplus
namespace CP2K_NS {
#endif
const int PLANEWAVES = 401;
const int GAUSSIAN = 402;
const int SQUARE = 391;
const int SQUAREROOT = 392;
const int PW_REALDATA1D = REALDATA1D;
const int PW_COMPLEXDATA1D = COMPLEXDATA1D;
const int PW_REALDATA3D = REALDATA3D;
const int PW_COMPLEXDATA3D = COMPLEXDATA3D;
typedef struct coeff_type {
struct pw_type pw;
int in_use;
} coeff_type;
void coeff_allocate_from_coeff(coeff_type, coeff_type, int);
void coeff_allocate_from_grid(coeff_type, pw_grid_type, int);
void coeff_deallocate(coeff_type);
void coeff_zero(coeff_type);
#ifdef __cplusplus
}
#endif
#endif

9
src/complex.h Normal file
View file

@ -0,0 +1,9 @@
#ifndef _COMPLEX_H
#define _COMPLEX_H
typedef struct complex {
double real;
double imag;
} complex;
#endif

934
src/constraint.F → src/constraint.c
View file

File diff suppressed because it is too large Load diff

48
src/convert_units.F → src/convert_units.c
View file

@ -1,6 +1,6 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
! Copyright (C) 2024 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE convert_units
@ -14,22 +14,22 @@ MODULE convert_units
USE physcon, ONLY : joule, evolt, boltzmann
USE simulation_cell, ONLY : cell_type
USE unit, ONLY : unit_convert_type
IMPLICIT NONE
PRIVATE
PUBLIC :: convert
CONTAINS
!-----------------------------------------------------------------------------!
!
! routine to convert input variables from external units to internal units
SUBROUTINE convert ( units, simpar, part, pstat, box, potparm, &
intra_param, ewald_param )
IMPLICIT NONE
! Arguments
TYPE (unit_convert_type ), INTENT ( IN ) :: units
TYPE (simulation_parameters_type ), INTENT ( INOUT ), OPTIONAL :: simpar
@ -39,27 +39,27 @@ SUBROUTINE convert ( units, simpar, part, pstat, box, potparm, &
TYPE (ewald_parameters_type ), INTENT ( INOUT ), OPTIONAL :: ewald_param
TYPE (potentialparm_type ), INTENT ( INOUT ), OPTIONAL :: potparm ( :, : )
TYPE (particle_prop_type ), INTENT ( INOUT ), OPTIONAL :: pstat ( : )
! Locals
REAL ( dbl ), PARAMETER :: permittivity = 8.854187817E-12_dbl
REAL ( dbl ), PARAMETER :: permittivity = 8.854187817E-12
REAL ( dbl ) :: conv2k
INTEGER :: i
!------------------------------------------------------------------------------
! convert hmat
IF ( PRESENT ( box)) box%hmat = box%hmat*units%lconv
IF ( PRESENT ( simpar)) THEN
! convert time
simpar%dt = simpar%dt*units%tconv
! convert temperature
simpar%temp_ext = simpar%temp_ext*units%econv
! convert pressure
! convert pressure
simpar%p_ext = simpar%p_ext/units%pconv
! convert verlet skin
simpar%verlet_skin = simpar%verlet_skin*units%lconv
! converting extended system time constants
SELECT CASE (simpar%ensemble)
CASE ( 'NVT')
@ -72,14 +72,14 @@ SUBROUTINE convert ( units, simpar, part, pstat, box, potparm, &
simpar%tau_nhc = simpar%tau_nhc*units%tconv
END SELECT
END IF
IF ( PRESENT ( pstat)) THEN
! convert mass
DO i = 1, size(pstat)
pstat(i) %mass = pstat(i) %mass*units%mconv
END DO
END IF
IF ( PRESENT ( part)) THEN
! convert position
DO i = 1, size(part)
@ -88,7 +88,7 @@ SUBROUTINE convert ( units, simpar, part, pstat, box, potparm, &
part(i) %r(3) = part(i) %r(3)*units%lconv
END DO
END IF
IF ( PRESENT ( potparm)) THEN
! convert LJ parameters
potparm ( :, : ) %lj%epsilon = potparm ( :, : ) %lj%epsilon*units%econv
@ -124,9 +124,9 @@ SUBROUTINE convert ( units, simpar, part, pstat, box, potparm, &
! tbhop(i,j) %gsp(:) %rt = tbhop(i,j) %gsp(:) %rt*units%lconv
! END DO
! END DO
!
! note: e_cutoff_coul and energy_cutoff do not need to be converted
! because ewald_correction and the initialization of the energy_cutoff
!
! note: e_cutoff_coul and energy_cutoff do not need to be converted
! because ewald_correction and the initialization of the energy_cutoff
! are called later and the cut-off will be computed in the appropriate units
!
! convert electrostatic parameters
@ -134,17 +134,17 @@ SUBROUTINE convert ( units, simpar, part, pstat, box, potparm, &
ewald_param % eps0 = permittivity*units%eps0
ewald_param % alpha = ewald_param % alpha/units % lconv
END IF
IF ( PRESENT ( intra_param)) THEN
! converting bond, bends and torsion and constraint parameters
intra_param % bond_param ( :, : ) % r0 = intra_param % bond_param ( :, : ) % r0 &
*units%lconv
intra_param % bond_param ( :, : ) % k = intra_param % bond_param ( :, : ) % k &
* units%econv/units%lconv/units%lconv
intra_param % bend_param(:,:,:) % k &
= intra_param % bend_param(:,:,:) % k * units%econv
DO i = 1, 6
intra_param % torsion_param(:,:,:,:) % parm(i) = &
intra_param % torsion_param(:,:,:,:) % parm(i)*units%econv
@ -152,7 +152,7 @@ SUBROUTINE convert ( units, simpar, part, pstat, box, potparm, &
intra_param % constraint_distance ( :, : ) &
= intra_param % constraint_distance ( :, : ) * units%lconv
END IF
END SUBROUTINE convert
!******************************************************************************

69
src/cp2k.F
View file

@ -1,69 +0,0 @@
!-----------------------------------------------------------------------------!
! !
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
! !
! This program is free software; you can redistribute it and/or modify !
! it under the terms of the GNU General Public License as published by !
! the Free Software Foundation; either version 2 of the License, or !
! (at your option) any later version. !
! !
! This program is distributed in the hope that it will be useful, !
! but WITHOUT ANY WARRANTY; without even the implied warranty of !
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !
! GNU General Public License for more details. !
! !
! You should have received a copy of the GNU General Public License !
! along with this program; if not, write to the Free Software !
! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. !
! !
!-----------------------------------------------------------------------------!
!
!! Main Program
!! - performs only basic initialisations
!! - calls major modules
PROGRAM cp2k
USE cp2k_input, ONLY : read_cp2k_section
USE environment, ONLY : initialisation, trailer
USE fist, ONLY : fist_main
USE fist_global, ONLY : set_fist_global
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl, print_kind_info
USE parallel, ONLY : start_parallel, end_parallel
USE physcon, ONLY : print_physcon
USE tbmd, ONLY : tbmd_main
USE tbmd_global, ONLY : set_tbmd_global
IMPLICIT NONE
TYPE ( global_environment_type ) :: globenv
!------------------------------------------------------------------------------
CALL start_parallel ( globenv )
CALL initialisation ( globenv )
CALL read_cp2k_section ( globenv )
IF ( globenv % ionode .AND. globenv % print_level > 9 ) &
CALL print_kind_info ( globenv % scr )
IF ( globenv % ionode .AND. globenv % print_level > 4 ) &
CALL print_physcon ( globenv % scr )
IF ( globenv % program_name == 'FIST' ) THEN
CALL set_fist_global ( globenv )
CALL fist_main()
ELSE IF ( globenv % program_name == 'TBMD' ) THEN
CALL set_tbmd_global ( globenv )
CALL tbmd_main()
END IF
CALL trailer ( globenv )
CALL end_parallel()
END PROGRAM cp2k

65
src/cp2k.cpp Normal file
View file

@ -0,0 +1,65 @@
/*---------------------------------------------------------------------------*/
/* */
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 2024 CP2K developers group */
/* */
/* This program is free software; you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation; either version 2 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program; if not, write to the Free Software */
/* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
/* */
/*---------------------------------------------------------------------------*/
#include "global_types.h"
// Main Program
// - performs only basic initialisations
// - calls major modules
using namespace CP2K_NS;
class cp2k_input;
class environment;
class fist;
class fist_global;
class parallel;
class physcon;
class tbmd;
class tbmd_global;
int cp2k(int argc, char **argv) {
global_environment_type globenv;
start_parallel(globenv);
initialisation(globenv);
read_cp2k_section(globenv);
if (globenv.ionode && globenv.print_level > 9) print_kind_info(globenv.scr);
if (globenv.ionode && globenv.print_level > 4) print_physcon(globenv.scr);
if (globenv.program_name == "FIST") {
set_fist_global(globenv);
fist_main();
} else if (globenv.program_name == "TBMD") {
set_tbmd_global(globenv);
tbmd_main();
}
trailer(globenv);
end_parallel();
}

14
src/cp2k.h Normal file
View file

@ -0,0 +1,14 @@
#ifndef _LAMMPS_H
#define _LAMMPS_H
#ifdef __cplusplus
namespace CP2K_NS {
#endif
int cp2k(int, char **);
#ifdef __cplusplus
}
#endif
#endif

114
src/cp2k_input.F
View file

@ -1,114 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE cp2k_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
USE string_utilities, ONLY : uppercase
PRIVATE
PUBLIC :: read_cp2k_section
CONTAINS
!!>---------------------------------------------------------------------------!
!! !
!! SECTION: &cp2k ... &end !
!! !
!! program prg_name !
!! project name !
!! iolevel n !
!! !
!!<---------------------------------------------------------------------------!
SUBROUTINE read_cp2k_section ( globenv )
IMPLICIT NONE
! Arguments
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
! Locals
INTEGER :: ierror, ilen, iw, source, allgrp
CHARACTER ( LEN = 6 ) :: string
CHARACTER ( LEN = 5 ) :: label
! ..defaults
globenv%program_name = 'NONE'
globenv%project_name = 'project'
globenv%print_level = 0
iw = globenv % scr
! ..parse the input section
label = '&CP2K'
ierror = 0
CALL parser_init(globenv%input_file_name,label,ierror,globenv)
IF (ierror /= 0) THEN
IF (globenv % ionode) &
WRITE ( iw, '( a )' ) ' No input section &CP2K found '
ELSE
CALL read_line()
DO WHILE ( test_next() /= 'X' )
ilen = 6
CALL cfield ( string, ilen )
CALL uppercase ( string )
SELECT CASE (string)
CASE DEFAULT
CALL p_error()
CALL stop_prg( 'input_cntl','unknown option')
CASE ( 'PROGRA' )
ilen = 0
CALL cfield ( globenv % program_name, ilen )
CALL uppercase ( globenv % program_name )
CASE ( 'PROJEC' )
ilen = 0
CALL cfield ( globenv % project_name, ilen )
CASE ( 'IOLEVE')
globenv%print_level = get_int()
END SELECT
CALL read_line()
END DO
END IF
CALL parser_end
! ..end of parsing the input section
IF ( globenv % program_name == 'QUICKSTEP' ) globenv % program_name = 'QS'
IF ( globenv % ionode ) THEN
! ..write some information to output
WRITE ( iw, '( A, T41, A )' ) ' CP2K| Program name ', &
ADJUSTR ( globenv % program_name ( 1:40 ) )
WRITE ( iw, '( A, T41, A )' ) ' CP2K| Input file name ', &
ADJUSTR ( globenv % input_file_name ( 1:40 ) )
WRITE ( iw, '( A, T41, A )' ) ' CP2K| Project name ', &
ADJUSTR ( globenv % project_name ( 1:40 ) )
WRITE ( iw, '( A, T77, I4 )' ) ' CP2K| Global print level ', &
globenv % print_level
WRITE ( iw, '( A, T75, I6 )' ) ' CP2K| Total number of processors ', &
globenv % num_pe
WRITE ( iw, '( A, T75, I6 )' ) ' CP2K| This output from processor ', &
globenv % mepos
WRITE (iw,*)
END IF
END SUBROUTINE read_cp2k_section
!******************************************************************************
END MODULE cp2k_input

89
src/cp2k_input.cpp Normal file
View file

@ -0,0 +1,89 @@
/*---------------------------------------------------------------------------*/
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 2024 CP2K developers group */
/*---------------------------------------------------------------------------*/
#include <string>
#include <iomanip>
#include <iostream>
#include "global_types.h"
using namespace CP2K_NS;
void read_cp2k_section(global_environment_type *globenv) {
// Locals
int ierror, ilen, iw, source, allgrp;
std::string string;
std::string label;
// defaults
globenv->program_name = "None";
globenv->project_name = "project";
globenv->print_level = 0;
iw = globenv->scr;
// parse the input section
label = "&CP2K";
ierror = 0;
parser_init(globenv->input_file_name, label, &ierror, globenv);
if (ierror != 0 ) {
if (globenv->ionode) std::cout << " No input section &CP2K found " << std::endl;
} else {
read_line()
while (test_next() != "X") {
ilen = 6;
cfield(string, ilen);
string.upper()
switch (string) {
case "PROGRA":
ilen = 0
cfield(globenv->program_name, ilen);
globenv->program_name.upper();
break;
case "PROJEC":
ilen = 0
cfield(globenv->project_name, ilen);
break;
case "IOLEVE":
globenv->print_level = get_int();
break;
default:
p_error();
stop_prg("input_cntl", "unknown option");
}
read_line();
}
}
parser_end();
// end of parsing the input section
if (globenv->program_name == "QUICKSTEP") globenv->program_name = "QS";
if (globenv->ionode) {
// write some information to output
std::cout << std::left << std::setw(40) << " CP2K| Program name ";
std::cout << std::right << std::setw(40) << globenv->program_name << std::endl;
std::cout << std::left << std::setw(40) << " CP2K| Input file name ";
std::cout << std::right << std::setw(40) << globenv->input_file_name << std::endl;
std::cout << std::left << std::setw(40) << " CP2K| Project name ";
std::cout << std::right << std::setw(40) << globenv->project_name << std::endl;
std::cout << std::left << std::setw(76) << " CP2K| Global print level ";
std::cout << globenv->print_level << std::endl;
std::cout << std::left << std::setw(74) << " CP2K| Total number of processors ";
std::cout << globenv->num_pe << std::endl;
std::cout << std::left << std::setw(74) << " CP2K| This output from processor ";
std::cout << globenv->mepos << "\n" << std::endl;
}
}

30
src/dg.h Normal file
View file

@ -0,0 +1,30 @@
#ifndef _DG_H
#define _DG_H
#include "coefficients.h"
#include "particle.h"
#ifdef __cplusplus
namespace CP2K_NS {
#endif
typedef struct dg_rho0_type {
char type[13];
int grid;
double *gcc;
double *zet;
coeff_type density;
} dg_rho0_type;
typedef struct dg_type {
int grid_index;
int nparts;
dg_rho0_type dg_rho0;
particle_list_type *plist;
} dg_type;
#ifdef __cplusplus
}
#endif
#endif

0
src/dg_rho0s.F → src/dg_rho0s.c
View file

0
src/dg_types.F → src/dg_types.c
View file

0
src/dgs.F → src/dgs.c
View file

0
src/dump.F → src/dump.c
View file

0
src/eigenvalueproblems.F → src/eigenvalueproblems.c
View file

176
src/environment.F
View file

@ -1,176 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE environment
USE global_types, ONLY : global_environment_type
USE kinds, ONLY : dbl
USE machine, ONLY : m_hostnm, m_getcwd, m_getlog, &
m_getuid, m_getpid, m_getarg
USE mathconstants, ONLY : zero
USE physcon, ONLY : init_physcon
USE stop_program, ONLY : set_stop, stop_prg
USE timesl, ONLY : walltime, cputime, datum
USE timings, ONLY : timeprint, timeset, timestop
USE util, ONLY : ran2
IMPLICIT NONE
PRIVATE
PUBLIC :: initialisation, trailer
CHARACTER ( LEN = 32 ) :: hname = ' ' !! name of the computer
CHARACTER ( LEN = 60 ) :: curdir = ' ' !! current directory name
CHARACTER ( LEN = 10 ) :: user = ' ' !! user name
INTEGER :: my_pid = 0 !! process ID
INTEGER :: my_uid = 0 !! user ID
INTEGER :: handle
CONTAINS
!******************************************************************************
SUBROUTINE initialisation ( globenv )
IMPLICIT NONE
! Arguments
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv ! check intent
! Locals
CHARACTER ( LEN = 26 ) :: datx
INTEGER :: iw, l, initial_random
!------------------------------------------------------------------------------
globenv % scr = 6
CALL set_stop ( globenv % scr )
IF ( globenv % ionode ) THEN
iw = globenv % scr
CALL set_environment ( globenv )
CALL datum ( datx )
WRITE ( iw, &
'( A, T55, A26, /, A, T49, A32, /, A, T71, A10, /, A, T71, I10 )' ) &
' **** **** ****** ** PROGRAM STARTED AT ', ADJUSTR ( datx ), &
' ***** ** *** *** ** PROGRAM STARTED ON ', ADJUSTR ( hname ), &
' ** **** ****** PROGRAM STARTED BY ', ADJUSTR ( user ), &
' ***** ** ** ** ** PROGRAM PROCESS ID ', my_pid
l = LEN_TRIM ( curdir )
IF ( l <= 38 ) THEN
WRITE ( iw, '( A, T43, A38 )' ) &
' **** ** ******* ** PROGRAM STARTED IN ', &
ADJUSTR ( curdir ( 1 : l ) )
ELSE
WRITE ( iw, '( A )' ) &
' **** ** ******* ** PROGRAM STARTED IN '
WRITE ( iw, '( T21, A60 )' ) ADJUSTR ( curdir )
END IF
WRITE ( iw, '( )' )
END IF
! initialize timing
CALL timeset ( 'CP2K', 'I', ' ', handle )
! initialize random number generator
globenv % idum = -1
initial_random = ran2 ( globenv % idum )
! initialize physical constants
CALL init_physcon()
END SUBROUTINE initialisation
!******************************************************************************
SUBROUTINE trailer ( globenv )
IMPLICIT NONE
! Arguments
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv ! check intent
! Locals
CHARACTER ( LEN = 26 ) :: datx
INTEGER :: iw, l
!------------------------------------------------------------------------------
CALL timestop ( zero, handle )
IF ( globenv % ionode ) THEN
iw = globenv % scr
CALL timeprint ( iw, globenv % print_level )
WRITE ( iw, '( )' )
CALL datum ( datx )
WRITE ( iw, &
'( A, T55, A26, /, A, T49, A32, /, A, T71, A10, /, A, T71, I10 )' ) &
' **** **** ****** ** PROGRAM ENDED AT ', ADJUSTR ( datx ), &
' ***** ** *** *** ** PROGRAM RAN ON ', ADJUSTR ( hname ), &
' ** **** ****** PROGRAM RAN BY ', ADJUSTR ( user ), &
' ***** ** ** ** ** PROGRAM PROCESS ID ', my_pid
l = LEN_TRIM ( curdir )
IF ( l <= 38 ) THEN
WRITE ( iw, '( A, T43, A38 )' ) &
' **** ** ******* ** PROGRAM RAN IN ', &
ADJUSTR ( curdir ( 1 : l ) )
ELSE
WRITE ( iw, '( A )' ) &
' **** ** ******* ** PROGRAM RAN IN '
WRITE ( iw, '( T21, A60 )' ) ADJUSTR ( curdir )
END IF
END IF
END SUBROUTINE trailer
!******************************************************************************
SUBROUTINE set_environment ( globenv )
IMPLICIT NONE
! Arguments
TYPE ( global_environment_type ), INTENT ( INOUT ) :: globenv
! Locals
INTEGER :: narg
INTEGER :: iargc
!------------------------------------------------------------------------------
CALL m_hostnm ( hname )
CALL m_getcwd ( curdir )
CALL m_getlog ( user )
CALL m_getuid ( my_uid )
CALL m_getpid ( my_pid )
narg = iargc()
SELECT CASE ( narg )
CASE ( 0 )
globenv % input_file_name = 'input'
CASE ( 1 )
CALL m_getarg ( 1, globenv % input_file_name )
CASE DEFAULT
CALL stop_prg ( 'set_environment', &
'illegal number of command line arguments passed' )
END SELECT
END SUBROUTINE set_environment
!******************************************************************************
END MODULE environment

120
src/environment.cpp Normal file
View file

@ -0,0 +1,120 @@
/*---------------------------------------------------------------------------*/
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 2024 CP2K developers group */
/*---------------------------------------------------------------------------*/
#include "global_types.h"
#ifdef __cplusplus
#include <string>
#include <iostream>
using namespace CP2K_NS;
#endif
void initialisation(global_environment_type *globenv) {
std::string datx;
int iw, l, initial_random;
//-----------------------------------------------------------------------------
globenv->scr = 6;
set_stop(globenv.scr);
if (globenv->ionode) {
iw = globenv.scr;
set_environment(globenv);
datum(datx);
std::cout << std::left << std::setw(54) << " **** **** ****** ** PROGRAM STARTED AT ";
std::cout << std::right << std::setw(26) << datx << std::endl;
std::cout << std::left << std::setw(48) << " ***** ** *** *** ** PROGRAM STARTED ON ";
std::cout << std::right << std::setw(32) << hname << std::endl;
std::cout << std::left << std::setw(70) << " ** **** ****** PROGRAM STARTED BY ";
std::cout << std::right << std::setw(10) << user << std::endl;
std::cout << std::left << std::setw(70) << " ***** ** ** ** ** PROGRAM PROCESS ID ";
std::cout << std::right << std::setw(10) << my_pid << std::endl;
std::cout << std::left << std::setw(42) << " **** ** ******* ** PROGRAM STARTED IN ";
std::cout << std::right << std::setw(38) << curdir << "\n" << std::endl;
}
// initialize timing
timeset("CP2K", "I", " ", handle);
// initialize random number generator
globenv.idum = -1;
initial_random = ran2(globenv.idum);
init_physcon();
}
void trailer(global_environment_type *globenv) {
std::string datx;
int iw, l;
timestop(zero, handle);
if (globenv.ionode) {
iw = globenv.scr;
timeprint(iw, globenv.print_level);
std::cout << "\n" << std::endl;
datum(datx);
std::cout << std::left << std::setw(54) << " **** **** ****** ** PROGRAM ENDED AT ";
std::cout << std::right << std::setw(26) << datx << std::endl;
std::cout << std::left << std::setw(48) << " ***** ** *** *** ** PROGRAM RAN ON ";
std::cout << std::right << std::setw(32) << hname << std::endl;
std::cout << std::left << std::setw(70) << " ** **** ****** PROGRAM RAN BY ";
std::cout << std::right << std::setw(10) << user << std::endl;
std::cout << std::left << std::setw(70) << " ***** ** ** ** ** PROGRAM PROCESS ID ";
std::cout << std::right << std::setw(10) << my_pid << std::endl;
std::cout << std::left << std::setw(70) << " **** ** ******* ** PROGRAM RAN IN ";
std::cout << std::right << std::setw(10) << curdir << std::endl;
}
}
/*****************************************************************************/
void set_environment(global_environment_type *globenv) {
// Locals
int narg, iargc;
//-----------------------------------------------------------------------------
m_hostnm(hname);
m_getcwd(curdir);
m_getlog(user);
m_getuid(my_uid);
m_getpid(my_pid);
narg = iargc();
switch (narg) {
case 0:
globenv.input_file_name = "input";
break;
case 1:
m_getarg(1, globenv.input_file_name);
break;
default:
stop_prg("set_environment", "illegal number of command line arguments passed")
}
}

0
src/erf.F → src/erf.c
View file

0
src/erf_fn.F → src/erf_fn.c
View file

23
src/ewald.h Normal file
View file

@ -0,0 +1,23 @@
#ifndef _EWALD_H
#define _EWALD_H
#ifdef __cplusplus
namespace CP2K_NS {
#endif
typedef struct ewald_parameters_type {
char ewald_type[13];
double alpha;
double eps0;
double epsilon;
int gmax;
int ns_max;
int gtot;
int ewald_grp;
} ewald_parameters_type;
#ifdef __cplusplus
}
#endif
#endif

0
src/ewald_parameters_types.F → src/ewald_parameters_types.c
View file

428
src/ewalds.F
View file

@ -1,428 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE ewalds
USE coefficient_types, ONLY : coeff_type
USE dgs, ONLY : dg_grid_setup
USE dg_types, ONLY : dg_type
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE kinds, ONLY : dbl
USE mathconstants, ONLY : pi, zero
USE md, ONLY : thermodynamic_type
USE molecule_types, ONLY : particle_node_type
USE mp, ONLY : mp_sum
USE particle_types, ONLY : particle_type
USE pme, ONLY : pme_setup
USE pw_grid_types, ONLY : pw_grid_type, HALFSPACE
USE pw_grids, ONLY : pw_find_cutoff, pw_grid_setup
USE simulation_cell, ONLY : cell_type
USE stop_program, ONLY : stop_prg, stop_memory
USE structure_factor_types, ONLY : structure_factor_type
USE structure_factors, ONLY : structure_factor_evaluate, &
structure_factor_allocate, &
structure_factor_deallocate
USE timings, ONLY : timeset, timestop
USE util, ONLY : matvec_3x3
IMPLICIT NONE
PRIVATE
PUBLIC :: ewald_gaussian, ewald_correction, ewald_self
PUBLIC :: ewald_print, ewald_initialize
INTEGER :: ewald_grp
CONTAINS
!******************************************************************************
#if defined ( __PGI )
#define mp mpp
#endif
!******************************************************************************
SUBROUTINE ewald_gaussian (dg, ewald_param, fg_coulomb, vg_coulomb, &
pv_g, pnode, box )
! computes the potential and the force from the g-space part of
! the 1/r potential
! Ref.: J.-P. Hansen, Enrico Fermi School, 1985
! Note: Only the positive G-vectors are used in the sum.
IMPLICIT NONE
! Arguments
TYPE ( dg_type ), INTENT ( IN ), TARGET :: dg
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
TYPE ( particle_node_type ), INTENT ( IN ), DIMENSION ( : ) :: pnode
REAL ( dbl ), INTENT ( OUT ) :: vg_coulomb
REAL ( dbl ), INTENT ( OUT ), DIMENSION ( :, : ) :: fg_coulomb
REAL ( dbl ), INTENT ( OUT ), DIMENSION ( :, : ) :: pv_g
! Locals
INTEGER :: node, gpt, nnodes, handle, isos, gptl, gptm, gptn
INTEGER :: lp, mp, np
INTEGER, DIMENSION ( :, : ), POINTER :: bds
REAL ( dbl ) :: gdotr, gsq, four_alpha_sq, common, pref
REAL ( dbl ) :: denom, charge, e_igdotr, flops
REAL ( dbl ) :: gauss
REAL ( dbl ), DIMENSION ( 3 ) :: vec
REAL ( dbl ), DIMENSION ( :, :, : ), POINTER :: rho0
COMPLEX ( dbl ), DIMENSION ( : ), ALLOCATABLE :: summe
TYPE ( coeff_type ), POINTER :: dg_rho0
TYPE ( pw_grid_type ), POINTER :: pw_grid
TYPE ( structure_factor_type ) :: exp_igr
!------------------------------------------------------------------------------
CALL timeset ( 'EWALD', 'E', 'Mflops', handle )
flops = 0.0_dbl
! pointing
dg_rho0 => dg % dg_rho0 % density
rho0 => dg % dg_rho0 % density % pw % cr3d
pw_grid => dg % dg_rho0 % density % pw % pw_grid
bds => dg % dg_rho0 % density % pw % pw_grid % bounds
! allocating
nnodes = SIZE ( pnode )
CALL structure_factor_allocate ( pw_grid % bounds, nnodes, exp_igr)
ALLOCATE ( summe ( 1:pw_grid % ngpts_cut ), STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'ewald', 'summe', 0 )
! Initializing vg_coulomb and fg_coulomb
vg_coulomb = 0.0_dbl
fg_coulomb = 0.0_dbl
pv_g = 0.0_dbl
! defining four_alpha_sq
four_alpha_sq = 4.0_dbl * ewald_param % alpha ** 2
DO node = 1, nnodes
vec = matvec_3x3 ( box % h_inv, pnode ( node ) % p % r )
CALL structure_factor_evaluate ( vec, pw_grid % npts, exp_igr % lb, &
exp_igr % ex(:,node), exp_igr % ey(:,node), exp_igr % ez(:,node) )
END DO
summe ( : ) = CMPLX ( 0.0_dbl, 0.0_dbl )
! looping over the positive g-vectors
DO gpt = 1, pw_grid % ngpts_cut
lp = pw_grid % mapl % pos ( pw_grid % g_hat ( 1, gpt ) )
mp = pw_grid % mapm % pos ( pw_grid % g_hat ( 2, gpt ) )
np = pw_grid % mapn % pos ( pw_grid % g_hat ( 3, gpt ) )
lp = lp + bds ( 1, 1 )
mp = mp + bds ( 1, 2 )
np = np + bds ( 1, 3 )
! initializing sum to be used in the energy and force
DO node = 1, nnodes
charge = pnode ( node ) % p % prop % charge
summe ( gpt ) = summe ( gpt ) + charge * &
( exp_igr % ex ( lp, node ) &
* exp_igr % ey ( mp, node ) &
* exp_igr % ez ( np, node ) )
END DO
END DO
flops = flops &
+ pw_grid % ngpts_cut * 10.0_dbl &
+ pw_grid % ngpts_cut * nnodes * 20.0_dbl
#if defined(__parallel)
CALL mp_sum ( summe, ewald_grp )
#endif
! looping over the positive g-vectors
DO gpt = 1, pw_grid % ngpts_cut
! computing the potential energy
lp = pw_grid % mapl % pos ( pw_grid % g_hat ( 1, gpt ) )
mp = pw_grid % mapm % pos ( pw_grid % g_hat ( 2, gpt ) )
np = pw_grid % mapn % pos ( pw_grid % g_hat ( 3, gpt ) )
lp = lp + bds ( 1, 1 )
mp = mp + bds ( 1, 2 )
np = np + bds ( 1, 3 )
IF ( pw_grid % gsq ( gpt ) <= 1.0E-10_dbl ) CYCLE
gauss = ( rho0 ( lp, mp, np ) * pw_grid % vol ) ** 2 &
/ pw_grid % gsq ( gpt )
common = gauss * REAL ( summe ( gpt ) * CONJG ( summe ( gpt ) ) )
vg_coulomb = vg_coulomb + common
! computing the force
DO node = 1, nnodes
e_igdotr = AIMAG ( summe(gpt) * CONJG &
( exp_igr % ex ( lp, node ) &
* exp_igr % ey ( mp, node ) &
* exp_igr % ez ( np, node ) ) )
charge = pnode ( node ) % p % prop % charge * gauss * e_igdotr
fg_coulomb ( :, node ) = fg_coulomb ( :, node ) &
+ charge * pw_grid % g ( :, gpt )
END DO
! compute the virial P*V
denom = 1.0_dbl / four_alpha_sq + 1.0_dbl / pw_grid % gsq ( gpt )
pv_g(1,1) = pv_g(1,1) + common * ( 1.0_dbl &
- 2.0_dbl * pw_grid % g ( 1, gpt ) * pw_grid % g ( 1, gpt ) * denom )
pv_g(1,2) = pv_g(1,2) - common * &
( 2.0_dbl * pw_grid % g ( 1, gpt ) * pw_grid % g ( 2, gpt ) * denom )
pv_g(1,3) = pv_g(1,3) - common * &
( 2.0_dbl * pw_grid % g ( 1, gpt ) * pw_grid % g ( 3, gpt ) * denom )
pv_g(2,1) = pv_g(2,1) - common * &
( 2.0_dbl * pw_grid % g ( 2, gpt ) * pw_grid % g ( 1, gpt ) * denom )
pv_g(2,2) = pv_g(2,2) + common * ( 1.0_dbl &
- 2.0_dbl * pw_grid % g ( 2, gpt ) * pw_grid % g ( 2, gpt ) * denom )
pv_g(2,3) = pv_g(2,3) - common * &
( 2.0_dbl * pw_grid % g ( 2, gpt ) * pw_grid % g ( 3, gpt ) * denom )
pv_g(3,1) = pv_g(3,1) - common * &
( 2.0_dbl * pw_grid % g ( 3, gpt ) * pw_grid % g ( 1, gpt ) * denom)
pv_g(3,2) = pv_g(3,2) - common * &
( 2.0_dbl * pw_grid % g ( 3, gpt ) * pw_grid % g ( 2, gpt ) * denom )
pv_g(3,3) = pv_g(3,3) + common * ( 1.0_dbl &
- 2.0_dbl * pw_grid % g ( 3, gpt ) * pw_grid % g ( 3, gpt ) * denom )
END DO
flops = flops + pw_grid % ngpts_cut * 55.0_dbl &
+ pw_grid % ngpts_cut * nnodes * 23.0_dbl
pref = 1.0_dbl / ( ewald_param % eps0 * box % deth )
vg_coulomb = vg_coulomb * pref
pv_g = pv_g * pref
pref = 2.0_dbl * pref
fg_coulomb = fg_coulomb * pref
CALL structure_factor_deallocate ( exp_igr )
DEALLOCATE ( summe, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'ewald', 'summe' )
flops = flops * 1.0E-6_dbl
CALL timestop ( flops, handle )
END SUBROUTINE ewald_gaussian
!******************************************************************************
!! computes the self interaction from g-space and the neutralizing background
SUBROUTINE ewald_self ( ewald_param, thermo, pnode )
! Arguments
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( IN ) :: pnode
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
! Locals
INTEGER :: nnodes, i
REAL ( dbl ) :: q_self, q_sum, q
!------------------------------------------------------------------------------
q_self = 0.0_dbl
q_sum = 0.0_dbl
nnodes = SIZE ( pnode )
DO i = 1, nnodes
q = pnode ( i ) % p % prop % charge
q_self = q_self + q*q
q_sum = q_sum + q
END DO
#if defined(__parallel)
CALL mp_sum(q_self,ewald_grp)
CALL mp_sum(q_sum,ewald_grp)
#endif
IF ( ewald_param % ewald_type == 'EWALD_GAUSS' .OR. &
ewald_param % ewald_type == 'PME_GAUSS' ) THEN
thermo % e_self = -q_self * ewald_param % alpha &
/ ( 4.0_dbl * ewald_param % eps0 * pi * SQRT ( pi ) )
thermo % e_neut = -q_sum ** 2 / ( 8.0_dbl * ewald_param % eps0 &
* ewald_param % alpha ** 2 )
END IF
END SUBROUTINE ewald_self
!******************************************************************************
!!> computing the shift in the potential energy
!!< for the real space coulomb potential
FUNCTION ewald_correction ( ewald_param, qi, qj, rcutsq ) &
RESULT ( e_cutoff_coul )
IMPLICIT NONE
! Arguments
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
REAL ( dbl ), INTENT ( IN ) :: qi, qj, rcutsq
REAL ( dbl ) ::e_cutoff_coul
! Locals
REAL ( dbl ) :: inv_fourpieps0
REAL ( dbl ), PARAMETER :: ac1 = 0.254829592_dbl, &
ac2 = -0.284496736_dbl, ac3 = 1.421413741_dbl, &
ac4 = -1.453152027_dbl, ac5 = 1.061405429_dbl, pc = 0.3275911_dbl
REAL ( dbl ) :: tc, arg, e_arg_arg, erfcf, erfc_over_rij
!------------------------------------------------------------------------------
! compute the value of the real-space coulomb at the cut-off
inv_fourpieps0 = 1.0_dbl / ( 4.0_dbl * pi * ewald_param % eps0 )
arg = ewald_param % alpha * SQRT ( rcutsq )
e_arg_arg = EXP ( - arg ** 2 )
tc = 1.0_dbl / ( 1.0_dbl + pc * arg )
erfcf = ((((ac5*tc+ac4)*tc+ac3)*tc+ac2)*tc+ac1)*tc * e_arg_arg
erfc_over_rij = erfcf / SQRT ( rcutsq )
e_cutoff_coul = inv_fourpieps0 * qi * qj * erfc_over_rij
END FUNCTION ewald_correction
!******************************************************************************
SUBROUTINE ewald_print ( iw, thermo, box, e_label )
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( IN ) :: iw
TYPE ( thermodynamic_type ), INTENT ( IN ) :: thermo
TYPE ( cell_type ), INTENT ( IN ) :: box
CHARACTER ( LEN = * ), INTENT ( IN ) :: e_label
!------------------------------------------------------------------------------
WRITE ( iw, '( A, A )' ) ' *********************************', &
'**********************************************'
WRITE ( iw, '( A, A, T35, A, T66, E15.7 )' ) ' INITIAL GSPACE ENERGY', &
e_label, '= ', thermo%gspace
WRITE ( iw, '( A, A, T35, A, T66, E15.7 )' ) ' SELF ENERGY CORRECTION', &
e_label, '= ', thermo%e_self
WRITE ( iw, '( A, A, T35, A, T66, E15.7 )' ) ' NEUT. BACKGROUND', e_label, &
'= ', thermo%e_neut/box%deth
WRITE ( iw, '( A, A, T35, A, T66, E15.7 )' ) ' BONDED CORRECTION', e_label, &
'= ', thermo%e_bonded
WRITE ( iw, '( A, A )' ) ' *********************************', &
'**********************************************'
END SUBROUTINE ewald_print
!******************************************************************************
SUBROUTINE ewald_initialize ( dg, part, pnode, pnode_grp, ewald_param, box, &
thermo, iounit, ewald_grid, pme_small_grid, pme_big_grid )
IMPLICIT NONE
! Arguments
TYPE ( dg_type ), INTENT ( OUT ) :: dg
TYPE ( ewald_parameters_type ), INTENT ( IN ) :: ewald_param
TYPE ( particle_type ), DIMENSION ( : ), INTENT ( IN ) :: part
TYPE ( particle_node_type ), DIMENSION ( : ), INTENT ( IN ) :: pnode
INTEGER, INTENT ( IN ) :: pnode_grp
TYPE ( cell_type ), INTENT ( IN ) :: box
TYPE ( thermodynamic_type ), INTENT ( INOUT ) :: thermo
INTEGER, INTENT ( IN ) :: iounit
TYPE ( pw_grid_type ), INTENT ( OUT ), OPTIONAL :: ewald_grid
TYPE ( pw_grid_type ), INTENT ( OUT ), OPTIONAL :: pme_small_grid
TYPE ( pw_grid_type ), INTENT ( OUT ), OPTIONAL :: pme_big_grid
! Locals
INTEGER :: natoms, iat, jat, gmax, npts_s ( 3 ), iw
REAL ( dbl ) :: qi, qj, i, cutoff
!------------------------------------------------------------------------------
! parallelisation is over atoms (pnodes), so the group of processors
! has to be the same as the group for the pnodes
ewald_grp = pnode_grp
! writing output to unit iounit
iw = iounit
natoms = SIZE ( part )
IF ( ewald_param % ewald_type /= 'NONE' ) THEN
WRITE ( iw, '( A,T71,A )' ) ' Ewald summation is done by:', &
ewald_param % ewald_type
WRITE ( iw, '( A,T71,F10.4 )' ) ' Ewald alpha parameter [A]', &
ewald_param % alpha
SELECT CASE ( ewald_param % ewald_type ( 1:3 ) )
CASE DEFAULT
WRITE ( iw, '( A,T71,I10 )' ) &
' Ewald G-space max. Miller index', ewald_param % gmax
CASE ( 'PME')
WRITE ( iw, '( A,T71,I10 )' ) &
' PME max small-grid points ', ewald_param % ns_max
WRITE ( iw, '( A,T71,F10.4 )' ) &
' PME gaussian tolerance ', ewald_param % epsilon
END SELECT
ELSE
WRITE ( iw, '( A )' ) ' No Ewald summation is performed'
END IF
! fire up the reciprocal space and compute self interaction and
! term from the neutralizing background.
IF ( ewald_param % ewald_type /= 'NONE' ) THEN
CALL ewald_self ( ewald_param, thermo, pnode )
! set up EWALD_GAUSS
IF ( PRESENT ( ewald_grid ) ) THEN
gmax = ewald_param % gmax
IF ( gmax == 2 * ( gmax / 2 ) ) THEN
CALL stop_prg ( "initialize_ewalds", "gmax has to be odd" )
END IF
ewald_grid % bounds ( 1, : ) = -gmax / 2
ewald_grid % bounds ( 2, : ) = +gmax / 2
npts_s = (/ gmax, gmax, gmax /)
ewald_grid % grid_span = HALFSPACE
CALL pw_find_cutoff ( npts_s, box, cutoff )
CALL pw_grid_setup( box, ewald_grid, cutoff)
CALL pme_setup (pnode, ewald_grid, ewald_param, dg )
END IF
! set up PME_GAUSS
IF ( PRESENT ( pme_small_grid ) .AND. PRESENT ( pme_big_grid ) ) THEN
npts_s ( : ) = ewald_param % ns_max
pme_small_grid % bounds ( 1, : ) = -npts_s ( : ) / 2
pme_small_grid % bounds ( 2, : ) = ( +npts_s ( : ) - 1 ) / 2
pme_small_grid % grid_span = HALFSPACE
pme_big_grid % grid_span = HALFSPACE
CALL dg_grid_setup ( box, npts_s, ewald_param % epsilon, &
ewald_param % alpha, pme_small_grid, &
pme_big_grid, ewald_param % ewald_type )
CALL pme_setup (pnode, pme_small_grid, ewald_param, dg )
END IF
END IF
END SUBROUTINE ewald_initialize
!******************************************************************************
END MODULE ewalds

330
src/ewalds.c Normal file
View file

@ -0,0 +1,330 @@
#include <complex.h>
#include <malloc.h>
#include <math.h>
#include <string.h>
#include "coefficients.h"
// #include "complex.h"
#include "dg.h"
#include "ewald.h"
#include "md.h"
#include "molecule.h"
#include "pw.h"
#include "simulation.h"
#include "structure.h"
#include "timings.h"
#include "util.h"
#if defined ( __PGI )
#define mp mpp
#endif
int ewald_grp;
void ewald_gaussian(dg_type dg, ewald_parameters_type ewald_param, double **fg_coulomb,
double vg_coulomb, double **pv_g, particle_node_type pnode[], cell_type box) {
//
int node, gpt, nnodes, handle, isos, gptl, gptm, gptn, lp, mp, np;
int *bds;
//
double gdotr, gsq, four_alpha_sq, common, pref, denom, charge, e_igdotr;
double flops, gauss;
size_t i, j;
double vec[3];
double ***rho0;
double complex *summe;
coeff_type *dg_rho0;
pw_grid_type *pw_grid;
structure_factor_type exp_igr;
timeset("EWALD", 'E', "Mflops", handle);
flops = 0.0;
// pointing
dg_rho0 = dg.dg_rho0->density;
rho0 = dg.dg_rho0->density->pw->cr3d;
pw_grid = dg.dg_rho0->density->pw->pw_grid;
bds = dg.dg_rho0.density.pw.pw_grid.bounds;
// allocating
nnodes = sizeof(pnode);
structure_factor_allocate(&pw_grid->bounds, nnodes, exp_igr);
summe = malloc(pw_grid->ngpts_cut*sizeof(double complex));
vg_coulomb = 0.0;
for (i = 0; i < node; i++) {
fg_coulomg[0][i] = 0.0;
fg_coulomg[1][i] = 0.0;
fg_coulomg[2][i] = 0.0;
}
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
pv_g[i][j] = 0.0;
}
}
// define four alpha_sq
four_alpha_sq = 4.0 * ewald_param.alpha * ewald_param.alpha;
for (node = 0; node < nnodes; node++) {
vec = matvec_3x3(box.h_inv, pnode[node].p->r);
structure_factor_evaluate(vec, pw_grid->npts, exp_igr.lb,
exp_igr.ex, exp_igr.ey, exp_igr.ez);
}
// looping over the positive g-vectors
for (gpt = 0; gpt < pw_grid->ngpts_cut; gpt++) {
lp = pw_grid->mapl->pos[pw_grid->g_hat[0][gpt]];
mp = pw_grid->mapm->pos[pw_grid->g_hat[1][gpt]];
np = pw_grid->mapn->pos[pw_grid->g_hat[2][gpt]];
lp += bds[0][0];
mp += bds[0][1];
np += bds[0][2];
for (node = 0; node < nnodes; node++) {
charge = pnode[node]->p->prop->charge;
summe[gpt] += charge *
( exp_igr.ex[lp][node]
* exp_igr.ey[mp][node]
* exp_igr.ez[np][node] );
}
}
flops += pw_grid->ngpts_cut * 10.0 + pw_grid->ngpts_cut * nnodes * 20.0;
#if defined (__parallel)
mp_sum(summe, ewald_grp);
#endif
// looping over the positive g-vectors
for (gpt = 0; gpt < pw_grid -> ngpts_cut; gpt++) {
// compute the potential energy
lp = pw_grid->mapl->pos[pw_grid->g_hat[0][gpt]];
mp = pw_grid->mapm->pos[pw_grid->g_hat[1][gpt]];
np = pw_grid->mapn->pos[pw_grid->g_hat[2][gpt]];
lp += bds[0][0];
mp += bds[0][1];
np += bds[0][2];
if (pw_grid->gsq[gpt] < 1.0E-10) break;
double tmp = rho0[lp][mp][np] * pw_grid->vol;
gauss = tmp * tmp / pw_grid->gsq[gpt];
common = gauss * (double)(summe[gpt] * conj(summe[gpt]));
vg_coulomb += common;
// computing the forces
for (node = 0; node < nnodes; node++) {
e_igdotr = cimag(summe[gpt] * conj
( exp_igr.ex[lp][node]
* exp_igr.ey[mp][node]
* exp_igr.ez[np][node]));
charge = pnode[node].p->prop->charge * gauss * e_igdotr;
fg_coulomb[0][node] += charge * pw_grid->g[0][gpt];
fg_coulomb[1][node] += charge * pw_grid->g[1][gpt];
fg_coulomb[2][node] += charge * pw_grid->g[2][gpt];
}
// computing the virial P*V
denom = 1.0 / four_alpha_sq + 1.0 / pw_grid->gsq[gpt];
pv_g[0][0] += common * (1.0 - 2.0 * pw_grid->g[0][gpt] * pw_grid->g[0][gpt] * denom);
pv_g[0][1] -= common * (2.0 * pw_grid->g[0][gpt] * pw_grid->g[1][gpt] * denom);
pv_g[0][2] -= common * (2.0 * pw_grid->g[0][gpt] * pw_grid->g[2][gpt] * denom);
pv_g[1][0] -= common * (2.0 * pw_grid->g[1][gpt] * pw_grid->g[0][gpt] * denom);
pv_g[1][1] += common * (1.0 - 2.0 * pw_grid->g[1][gpt] * pw_grid->g[1][gpt] * denom);
pv_g[1][2] -= common * (2.0 * pw_grid->g[1][gpt] * pw_grid->g[2][gpt] * denom);
pv_g[2][0] -= common * (2.0 * pw_grid->g[2][gpt] * pw_grid->g[0][gpt] * denom);
pv_g[2][1] -= common * (2.0 * pw_grid->g[2][gpt] * pw_grid->g[1][gpt] * denom);
pv_g[2][2] += common * (1.0 - 2.0 * pw_grid->g[2][gpt] * pw_grid->g[2][gpt] * denom);
}
flops += pw_grid->ngpts_cut * 55.0 + pw_grid->ngpts_cut * nnodes * 23.0;
pref = 1.0 / (ewald_param.eps0 * box.deth);
vg_coulomb *= pref;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
pv_g[i][j] *= pref;
}
}
pref *= 2.0;
for (j = 0; j < 3; j++) {
for (i = 0; i < nnodes; i++) {
fg_coulomb[j][i] *= pref;
}
}
for (size_t i = 0; i < nnodes; i++) {
fg_coulomb[0][i] *= pref;
fg_coulomb[1][i] *= pref;
fg_coulomb[2][i] *= pref;
}
structure_factor_deallocate(exp_igr);
free(summe);
flops *= 1.0E-6;
timestop(flops, handle);
}
void ewald_self(ewald_parameters_type ewald_param, thermodynamic_type *thermo, partilce_node_type pnode[], int nnodes) {
size_t i;
double q_self, q_sum, q;
q_self = 0.0;
q_sum = 0.0;
for (i = 0; i < nnodes; i++) {
q = pnode[i].p->prop->charge;
q_self += q * q;
q_sum += q;
}
#if defined (__parallel)
mp_sum(q_self, ewald_grp);
mp_sum(q_sum, ewald_grp);
#endif
if (strncmp(ewald_param.ewald_type, "EWALD_GAUSS", 20) ||
strncmp(ewald_param.ewald_type, "PME_GAUSS", 20)) {
thermo->e_self = -q_self * ewald_param.alpha
/ (4.0 * ewald_param.eps0 * sqrt(M_PI) * M_PI);
thermo->e_neut = -q_sum * q_sum / (8.0 * ewald_param.eps0
* ewald_param.alpha * ewald_param.alpha);
}
}
double ewald_correction(ewald_parameters_type ewald_param, double qi, double qj, double rcutsq) {
double e_cutoff_coul, inv_fourpiesp0;
double tc, arg, e_arg_arg, erfcf, erfc_over_rij;
double ac1 = 0.254829592;
double ac2 = -0.284496736;
double ac3 = 1.421413741;
double ac4 = -1.453152027;
double ac5 = 1.061405429;
double pc = 0.3275911;
// compute the value of the real-space coulomb at the cut-off
inv_fourpiesp0 = 1.0 / (4.0 * M_PI * ewald_param.eps0);
arg = ewald_param.alpha * sqrt(rcutsq);
e_arg_arg = exp(-arg * arg);
tc = 1.0 / (1.0 + pc * arg);
erfcf = ((((ac5 * tc + ac4) * tc + ac3) * tc + ac2) * tc + ac1) * tc * e_arg_arg;
erfc_over_rij = erfcf / sqrt(rcutsq);
e_cutoff_coul = inv_fourpiesp0 * qi * qj * erfc_over_rij;
return e_cutoff_coul;
}
void ewald_print(FILE *iw, thermodynamic_type thermo, cell_type box, char *e_label) {
fprintf(iw, "( %s, %s )", " *********************************",
"**********************************************");
fprintf(iw, "( %s, %s, %35s, %s, %66s, %15.7E )", " INITIAL GSPACE ENERGY",
e_label, "= ", thermo.gspace);
fprintf(iw, "( %s, %s, %35s, %s, %66s, %15.7E )", " SELF ENERGY CORRECTION",
e_label, "= ", thermo.e_self);
fprintf(iw, "( %s, %s, %35s, %s, %66s, %15.7E )", " NEUT. BACKGROUND", e_label,
"= ", thermo.e_neut / box.deth);
fprintf(iw, "( %s, %s, %35s, %s, %66s, %15.7E )", " BONDED CORRECTION", e_label,
"= ", thermo.e_bonded);
fprintf(iw, "( %s, %s )", " *********************************",
"**********************************************");
}
void ewald_initialize(dg_type *dg, particle_type *part, particle_node_type *pnode,
int pnode_grp, ewald_parameters_type *ewald_param, cell_type *box,
thermodynamic_type *thermo, int iounit, pw_grid_type *ewald_grid,
pw_grid_type *pme_small_grid, pw_grid_type *pme_big_grid) {
int natoms, iat, jat, gmax, npts_s[3], iw;
double qi, qj, i, cutoff;
// parallelisation is over atoms (pnodes), so the group of processors
// has to be the same as the group for the pnodes
ewald_grp = pnode_grp;
// writing output to unit iounit
iw = iounit;
natoms = sizeof(part) / sizeof(part[0]);
if (strcmp(ewald_param.ewald_type, "NONE") != 0) {
printf("Ewald summation is done by: %s\n", ewald_param.ewald_type);
printf("Ewald alpha parameter [A]: %f\n", ewald_param.alpha);
if (strncmp(ewald_param.ewald_type, "EWALD", 5) == 0) {
printf("Ewald G-space max. Miller index: %d\n", ewald_param.gmax);
} else if (strncmp(ewald_param.ewald_type, "PME", 3) == 0) {
printf("PME max small-grid points: %d\n", ewald_param.ns_max);
}
} else {
printf("No Ewald summation is performed\n");
}
// fire up the reciprocal space and compute self interaction and
// term from the neutralizing background.
if (strcmp(ewald_param.ewald_type, "NONE") != 0) {
ewald_self(ewald_param, thermo, pnode);
// set up EWALD_GAUSS
if (ewald_grid != NULL) {
gmax = ewald_param.gmax;
if (gmax % 2 == 0) {
stop_prg("initialize_ewalds", "gmax has to be odd");
}
ewald_grid->bounds[0][:] = -gmax / 2;
ewald_grid->bounds[1][:] = +gmax / 2;
npts_s[0] = gmax;
npts_s[1] = gmax;
npts_s[2] = gmax;
ewald_grid->grid_span = HALFSPACE;
pw_find_cutoff(npts_s, box, cutoff);
pw_grid_setup(box, ewald_grid, cutoff);
pme_setup(pnode, ewald_grid, ewald_param, dg);
}
// set up PME_GAUSS
if (pme_small_grid != NULL and pme_big_grid != NULL) {
npts_s[:] = ewald_param.ns_max;
pme_small_grid->bounds[0][:] = -npts_s[:]/2;
pme_small_grid->bounds[1][:] = (+npts_s[:] - 1) / 2;
pme_small_grid->grid_span = HALFSPACE;
pme_big_grid->grid_span = HALFSPACE;
dg_grid_setup(box, npts_s, ewald_param.epsilon, ewald_param.alpha,
pme_small_grid, pme_big_grid, ewald_param.ewald_type);
pme_setup(pnode, pme_small_grid, ewald_param, dg);
}
}

0
src/fft_tools.F → src/fft_tools.c
View file

328
src/fist.F
View file

@ -1,328 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE fist
USE atoms_input, ONLY : read_coord_vel, system_type
USE convert_units, ONLY : convert
USE dump, ONLY : dump_variables
USE ewalds, ONLY : ewald_print, ewald_correction
USE ewald_parameters_types, ONLY : ewald_parameters_type
USE header, ONLY : fist_header
USE kinds, ONLY : dbl
USE fist_debug, ONLY : debug_control
USE fist_input, ONLY : read_fist_section
USE fist_global, ONLY : fistpar
USE force_fields, ONLY : read_force_field_section, ATOMNAMESLENGTH
USE initialize_extended_types, ONLY : initialize_extended_type
USE initialize_molecule_types, ONLY : initialize_molecule_type
USE initialize_particle_types, ONLY : initialize_particle_type
USE integrator, ONLY : velocity_verlet, force, set_energy_parm, energy, &
set_integrator
USE input_types, ONLY : setup_parameters_type
USE linklist_control, ONLY : set_ll_parm
USE mathconstants, ONLY : zero
USE md, ONLY : read_md_section, simulation_parameters_type, &
initialize_velocities, thermodynamic_type, mdio_parameters_type
USE molecule_input, ONLY : read_molecule_section, read_setup_section
USE molecule_types, ONLY : molecule_type, intra_parameters_type
USE nose, ONLY : extended_parameters_type
USE pair_potential, ONLY : potentialparm_type, &
potential_f, spline_nonbond_control
USE particle_types, ONLY : particle_prop_type, particle_type
USE simulation_cell, ONLY : cell_type, get_hinv
USE stop_program, ONLY : stop_prg, stop_memory
USE structure_types, ONLY : structure_type
USE timings, ONLY : timeset, timestop
USE unit, ONLY : unit_convert_type, set_units
USE util, ONLY : close_unit, get_share
IMPLICIT NONE
PRIVATE
PUBLIC :: fist_main
TYPE ( mdio_parameters_type ) :: mdio
CONTAINS
!-----------------------------------------------------------------------------!
! FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST !
!-----------------------------------------------------------------------------!
SUBROUTINE fist_main
IMPLICIT NONE
! Locals
REAL ( dbl ) :: cons, ecut, qi, qj
REAL ( dbl ), DIMENSION ( :, : ), ALLOCATABLE :: rcut
INTEGER :: handle, itimes, isos, para_comm_cart
INTEGER :: iat, jat
CHARACTER ( LEN = ATOMNAMESLENGTH ), DIMENSION ( : ), POINTER :: atom_names
CHARACTER ( LEN = 20 ) :: set_fn
TYPE ( potentialparm_type ), DIMENSION ( :, : ), POINTER :: potparm
TYPE ( particle_prop_type ), DIMENSION ( : ), POINTER :: pstat
TYPE ( molecule_type ), DIMENSION ( : ), POINTER :: mol_setup
TYPE ( unit_convert_type ) :: units
TYPE ( simulation_parameters_type ) :: simpar
TYPE ( structure_type ) :: struc
TYPE ( extended_parameters_type ) :: nhcp
TYPE ( thermodynamic_type ) :: thermo
TYPE ( system_type ) :: ainp
TYPE ( ewald_parameters_type ) :: ewald_param
TYPE ( setup_parameters_type ) :: setup
TYPE ( intra_parameters_type ) :: intra_param
!------------------------------------------------------------------------------
CALL timeset ( 'FIST', 'I', ' ', handle )
IF ( fistpar % ionode ) CALL fist_header ( fistpar % scr )
!..read fist section
CALL read_fist_section ( setup, ewald_param, fistpar )
! read from the setup and molecule section of *.set
set_fn = setup % set_file_name
CALL read_setup_section ( mol_setup, set_fn, fistpar )
CALL read_molecule_section ( mol_setup, set_fn, fistpar )
! read force_field information for classical MD
CALL read_force_field_section ( setup, mol_setup, set_fn, &
intra_param, potparm, atom_names, pstat, fistpar )
!..read the input of the molecular dynamics section
CALL read_md_section ( simpar, fistpar, mdio )
simpar % program = 'FIST'
!..read atomic coordinates, velocities (optional) and the simulation box
ainp % rtype = simpar % read_type
CALL read_coord_vel ( ainp, setup % input_file_name, fistpar )
!..initialize box, perd
struc % box % hmat = ainp % box
struc % box % perd = setup % perd
CALL get_hinv ( struc % box )
!..initialize working units
CALL set_units ( setup % unit_type, units )
CALL set_energy_parm ( units % pconv, units % econv, units % l_label, &
units % vol_label, units % e_label, units % pres_label, &
units % temp_label, units % angl_label )
!..allocate memory for atoms and molecules
CALL allocmem ( ainp, mol_setup, struc )
!..initialize particle_type
CALL initialize_particle_type ( atom_names, simpar, mol_setup, &
ainp, pstat, struc % part )
!..convert the units
CALL convert ( units = units, simpar = simpar, &
part = struc % part, pstat = pstat, box = struc % box, &
potparm = potparm, intra_param = intra_param, &
ewald_param = ewald_param )
!..initialize molecule_type
CALL initialize_molecule_type ( mol_setup, intra_param, struc % pnode, &
struc % part, struc % molecule, fistpar )
!..initialize extended_parameters_type and get number of degrees of freedom
CALL initialize_extended_type ( struc % box, simpar, &
struc % molecule, mol_setup, nhcp, fistpar )
! initialize velocities if needed
IF ( simpar % read_type == 'POS' ) THEN
CALL initialize_velocities ( simpar, struc % part, fistpar )
END IF
!...initialize splines
potparm ( :, : ) % energy_cutoff = 0.0_dbl
potparm ( :, : ) % e_cutoff_coul = 0.0_dbl
CALL spline_nonbond_control ( potparm, pstat, 5000, ewald_param )
!..set linklist control parameters
ALLOCATE ( rcut ( setup % natom_type, setup % natom_type ), STAT = isos )
IF ( isos /=0 ) CALL stop_memory ( 'fist', 'rcut', 0 )
rcut ( :, : ) = potparm ( :, : ) % rcutsq
CALL set_ll_parm ( fistpar, simpar % verlet_skin, &
setup % natom_type, rcut, simpar % n_cell )
CALL set_ll_parm ( fistpar, printlevel = fistpar % print_level, &
ltype = 'NONBOND' )
DEALLOCATE ( rcut, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'fist', 'rcut' )
!..deallocate arrays needed for atom input
IF ( ASSOCIATED ( ainp % c ) ) THEN
DEALLOCATE ( ainp % c, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'fist', 'ainp%c' )
END IF
IF ( ASSOCIATED ( ainp % v ) ) THEN
DEALLOCATE ( ainp % v, STAT = isos )
IF ( isos /= 0 ) CALL stop_memory ( 'fist', 'ainp%v' )
END IF
!..initialize integrator
CALL set_integrator (fistpar % ionode, fistpar % group, fistpar % scr, mdio )
IF ( setup % run_type == 'DEBUG' ) THEN
!..debug the forces
CALL debug_control( fistpar, ewald_param, struc%part, struc%pnode, &
struc%molecule, struc%box, thermo, potparm, simpar % ensemble )
ELSE
!..MD
itimes = 0
CALL force ( struc, potparm, thermo, simpar, ewald_param )
IF ( fistpar % ionode .AND. ewald_param % ewald_type /= 'NONE' ) &
CALL ewald_print ( fistpar % scr, thermo, struc % box, &
units % e_label )
CALL energy ( itimes, cons, simpar, struc, thermo, nhcp )
DO itimes = 1, simpar % nsteps
CALL velocity_verlet ( itimes, cons, simpar, potparm, thermo, &
struc, ewald_param, nhcp )
IF ( MOD ( itimes, mdio % idump ) == 0 ) &
CALL dump_variables ( struc, mdio % dump_file_name, fistpar )
END DO
CALL dump_variables ( struc, mdio % dump_file_name, fistpar )
IF ( fistpar % ionode ) CALL close_unit ( 10, 99 )
END IF
!..deallocate memory for atoms and molecules
CALL deallocmem ( struc )
CALL timestop ( zero, handle )
END SUBROUTINE fist_main
!-----------------------------------------------------------------------------!
! FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST !
!-----------------------------------------------------------------------------!
SUBROUTINE allocmem ( ainp, mol_setup, struc )
IMPLICIT NONE
! Arguments
TYPE ( molecule_type ), DIMENSION ( : ), INTENT ( IN ) :: mol_setup
TYPE ( structure_type ), INTENT ( INOUT ) :: struc
TYPE ( system_type ), INTENT ( IN ) :: ainp
! Locals
INTEGER :: iw, natoms, nnodes, nmol, nmoltype, ios, iat, i
!------------------------------------------------------------------------------
struc % name = 'FIST MOLECULAR SYSTEM'
IF ( fistpar % num_pe == 1 ) THEN
natoms = SIZE ( ainp % c ( 1, : ) )
ALLOCATE ( struc % part ( 1:natoms ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'part', natoms )
ALLOCATE ( struc % pnode ( 1:natoms ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'pnode', natoms )
nmol = SUM ( mol_setup ( : ) % num_mol )
ALLOCATE ( struc % molecule ( 1:nmol ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'molecule', nmol )
IF ( fistpar % ionode .AND. fistpar % print_level > 3 ) THEN
iw = fistpar % scr
WRITE ( iw, '( A )' )
WRITE ( iw, '( A, T71, I10 )' ) &
' FIST| Number of allocated particles ', natoms
WRITE ( iw,'( A, T71, I10 )' ) &
' FIST| Number of allocated particle nodes ', natoms
WRITE ( iw, '( A, T71, I10 )' ) &
' FIST| Number of allocated molecules ', nmol
WRITE ( iw, '( A )' )
END IF
ELSE
!..replicated data
natoms = SIZE ( ainp % c ( 1, : ) )
ALLOCATE ( struc % part ( 1:natoms ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'part', natoms )
nmoltype = SIZE ( mol_setup )
nmol = 0
nnodes = 0
DO i = 1, nmoltype
nmol = nmol + get_share ( mol_setup ( i ) % num_mol, &
fistpar % num_pe, fistpar % mepos )
nnodes = nnodes + nmol * mol_setup ( i ) % molpar % natom
END DO
ALLOCATE ( struc % molecule ( 1:nmol ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'molecule' , nmol )
ALLOCATE ( struc % pnode ( 1:nnodes ), STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'pnode', nnodes )
IF ( fistpar % ionode .AND. fistpar % print_level > 3 ) THEN
iw = fistpar % scr
WRITE ( iw, '( A )' )
WRITE ( iw, '( A, T71, I10 )' ) &
' FIST| Number of allocated particles ', natoms
WRITE ( iw, '( A, T71, I10 )' ) &
' FIST| Number of allocated particle nodes ', nnodes
WRITE ( iw, '( A, I5, T71, I10 )' ) &
' FIST| Number of allocated molecules on processor ', &
fistpar % mepos, nmol
WRITE ( iw, '( A )' )
END IF
END IF
END SUBROUTINE allocmem
!******************************************************************************
SUBROUTINE deallocmem ( struc )
IMPLICIT NONE
! Arguments
TYPE ( structure_type ), INTENT ( INOUT ) :: struc
! Locals
INTEGER :: ios
!------------------------------------------------------------------------------
DEALLOCATE ( struc % part, STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'part' )
DEALLOCATE ( struc % pnode, STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'pnode' )
DEALLOCATE ( struc % molecule, STAT = ios )
IF ( ios /= 0 ) CALL stop_memory ( 'fist', 'molecule' )
END SUBROUTINE deallocmem
!******************************************************************************
END MODULE fist

241
src/fist.cpp Normal file
View file

@ -0,0 +1,241 @@
/*---------------------------------------------------------------------------*/
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 2024 CP2K developers group */
/*---------------------------------------------------------------------------*/
#include <vector>
#include "atoms_input.h"
#include "convert_units.h"
#include "dump.h"
#include "ewalds.h"
#include "ewald_parameters_types.h"
#include "header.h"
#include "kinds.h"
#include "fist_debug.h"
#include "fist_input.h"
#include "fist_global.h"
#include "force_fields.h"
#include "initialize_extended_types.h"
#include "initialize_molecule_types.h"
#include "initialize_particle_types.h"
#include "integrator.h"
#include "input_types.h"
#include "linklist_control.h"
#include "mathconstants.h"
#include "md.h"
#include "molecule_input.h"
#include "molecule_types.h"
#include "nose.h"
#include "pair_potential.h"
#include "particle_types.h"
#include "simulation_cell.h"
#include "stop_program.h"
#include "structure_types.h"
#include "timings.h"
#include "unit.h"
#include "util.h"
using namespace CP2K_NS;
class fist {
public:
void fist_main();
private:
void allocmem(system_type ainp, std::vector<molecule_type> mol_setup, structure_type *struc);
void deallocmem(structure_type struc);
};
/*----------------------------------------------------------------------------*/
/* FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST */
/*----------------------------------------------------------------------------*/
void fist::fist_main() {
double cons, ecut, qi, qj;
std::vector<std::vector<double>> rcut;
int handle, itimes, isos, para_comm_cart;
int iat, jat;
std::vector<std::string> atom_names;
std::string set_fn;
potentialparm_type potparm;
particle_prop_type pstat;
molecule_type mol_setup;
unit_convert_type units;
simulation_parameters_type simpar;
structure_type struc;
extended_parameters_type nhcp;
thermodynamic_type thermo;
system_type ainp;
ewald_parameters_type ewald_param;
setup_parameters_type setup;
intra_parameters_type intra_param;
timeset("FIST", "I", " ", handle);
if (fistpar.ionone) fist_header(fistpar.scr);
// read fist section
read_fist_section(setup, ewald_param, fistpar);
// read from the setup and molecule section of *.set
set_fn = setup.set_file_name;
read_setup_section(mol_setup, set_fn, fistpar);
read_molecule_section(mol_setup, set_fn, fistpar);
// read force_field information for classical MD
read_force_field_section(setup, mol_setup, set_fn,
intra_param, potparm, atom_names, pstat, fistpar);
// read the input of the molecular dynamics section
read_md_section(simpar, fistpar, mdio);
simpar.program = "FIST";
// read atomic coordinates, velocities (optional) and the simulation box
ainp.rtype = simpar.read_type;
atoms_input.read_coord_vel(ainp, setup.input_file_name, fistpar);
// initialize box, perd
struc.box.hmat = ainp.box;
struc.box.perd = setup.perd;
get_hinv(struc.box);
// initialize working units
set_units ( setup.unit_type, units );
set_energy_parm ( units.pconv, units.econv, units.l_label,
units.vol_label, units.e_label, units.pres_label,
units.temp_label, units.angl_label );
// allocate memory for atoms and molecules
allocmem ( ainp, mol_setup, struc );
// initialize particle_type
initialize_particle_type ( atom_names, simpar, mol_setup,
ainp, pstat, struc.part );
// convert the units
convert ( units = units, simpar = simpar,
part = struc.part, pstat = pstat, box = struc.box,
potparm = potparm, intra_param = intra_param,
ewald_param = ewald_param );
// initialize molecule_type
initialize_molecule_type ( mol_setup, intra_param, struc.pnode,
struc.part, struc.molecule, fistpar );
// initialize extended_parameters_type and get number of degrees of freedom
initialize_extended_type ( struc.box, simpar,
struc.molecule, mol_setup, nhcp, fistpar );
// initialize velocities if needed
if ( simpar.read_type == "POS" ) {
initialize_velocities ( simpar, struc.part, fistpar );
}
// initialize splines
potparm ( :, : ).energy_cutoff = 0.0;
potparm ( :, : ).e_cutoff_coul = 0.0;
spline_nonbond_control ( potparm, pstat, 5000, ewald_param );
// set linklist control parameters
allocate ( rcut ( setup.natom_type, setup.natom_type ), STAT = isos );
if ( isos !=0 ) stop_memory ( "fist", "rcut", 0 );
rcut ( :, : ) = potparm ( :, : ).rcutsq;
set_ll_parm ( fistpar, simpar.verlet_skin,
setup.natom_type, rcut, simpar.n_cell );
set_ll_parm ( fistpar, printlevel = fistpar.print_level,
ltype = "NONBOND" );
deallocate ( rcut, STAT = isos );
if ( isos != 0 ) stop_memory ( "fist", "rcut" );
// deallocate arrays needed for atom input
if (associated(ainp.c)) {
deallocate(ainp.c, STAT=isos);
if (isos != 0) stop_memory("fist", "ainp.c");
}
if (associated(ainp.v)) {
deallocate(ainp.v, STAT=isos);
if (isos != 0) stop_memory("fist", "ainp.v");
}
// initialize integrator
set_integrator(fistpar.ionode, fistpar.group, fistpar.scr, mdio);
if (setup.run_type == "DEBUG") {
// debug the forces
debug_control( fistpar, ewald_param, struc.part, struc.pnode,
struc.molecule, struc.box, thermo, potparm, simpar.ensemble );
} else {
// MD
itimes = 0;
force(struc, potparm, thermo, simpar, ewald_param);
if (fistpar.ionode && ewald_param.ewald_type != "None") {
ewald_print(fistpar.scr, thermo, struc.box, units.e_label);
}
energy(itimes, cons, simpar, struc, thermo, nhcp);
for (itimes = 0; itimes < simpar.nsteps; itimes++) {
velocity_verlet(itimes, cons, simpar, potparm, thermo,
struc, ewald_param, nhcp);
if (itimes % mdio.idump == 0) dump_variables(struc, mdio.dump_file_name, fistpar);
}
dump_variables(struc, mdio.dump_file_name, fistpar);
if (fistpar.ionode) close_unit(10,99);
}
// deallocate memory for atoms and molecules
deallocmem(struc);
timestop(zero, handle);
}
fist::~fist() {
deallocmem(struc);
}
/*----------------------------------------------------------------------------*/
/* FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST FIST */
/*----------------------------------------------------------------------------*/
void fist::allocmem(system_type ainp, std::vector<molecule_type> mol_setup, structure_type *struc) {
// Locals
int iw, natoms, nnodes, nmol, nmoltype, ios, iat, i;
//-----------------------------------------------------------------------------
struc.name = "FIST MOLECULAR SYSTEM"
if (fistpar.num_pe == 1) {
natoms = ainp.c.length();
allocate(struc.part, STAT = ios);
}
}
void fist::deallocmem(structure_type *struc) {
free(struc->part);
free(struc->pnode);
free(struc->molecule);
}

0
src/fist_debug.F → src/fist_debug.c
View file

0
src/fist_force.F → src/fist_force.c
View file

0
src/fist_force_numer.F → src/fist_force_numer.c
View file

0
src/fist_global.F → src/fist_global.c
View file

0
src/fist_input.F → src/fist_input.c
View file

0
src/fist_intra_force.F → src/fist_intra_force.c
View file

0
src/fist_nonbond_force.F → src/fist_nonbond_force.c
View file

0
src/force_control.F → src/force_control.c
View file

0
src/force_fields.F → src/force_fields.c
View file

0
src/force_fit.F → src/force_fit.c
View file

29
src/global_types.F
View file

@ -1,29 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE global_types
USE kinds, ONLY : dbl
IMPLICIT NONE
PRIVATE
PUBLIC :: global_environment_type
TYPE global_environment_type
CHARACTER ( LEN = 200 ) :: input_file_name
CHARACTER ( LEN = 200 ) :: program_name
CHARACTER ( LEN = 200 ) :: project_name
LOGICAL :: ionode !! this processor performs global output
INTEGER :: group !! group handle for all processors
INTEGER :: source !! processor id for global output node
INTEGER :: num_pe !! total number of processors
INTEGER :: mepos !! position number for local processor
INTEGER :: print_level !! how much information is written
INTEGER :: idum !! random number seed
INTEGER :: scr !! output unit
END TYPE global_environment_type
END MODULE global_types

32
src/global_types.h Normal file
View file

@ -0,0 +1,32 @@
#ifndef _GLOBAL_TYPES_H
#define _GLOBAL_TYPES_H
/*---------------------------------------------------------------------------*/
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 2024 CP2K developers group */
/*---------------------------------------------------------------------------*/
#ifdef __cplusplus
namespace CP2K_NS {
#include <string>
#endif
typedef struct global_environment_type {
char *input_file_name;
char *program_name;
char *project_name;
bool ionode; // this processor performs global output
int group; // group handle for all processors
int source; // processor id for global output node
int num_pe; // total number of processors
int mepos; // position number for local processor
int print_level; // how much information is written
int idum; // random number seed
int scr; // output unit
} global_environment_type;
#ifdef __cplusplus
}
#endif
#endif

10
src/header.c
View file

@ -1,6 +1,6 @@
/*---------------------------------------------------------------------------*/
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 2000 CP2K developers group */
/* Copyright (C) 2024 CP2K developers group */
/*---------------------------------------------------------------------------*/
#include "header.h"
@ -20,7 +20,7 @@ void fist_header(FILE *iw) {
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, " S. Balasubramanian / | 1998-2024\n");
fprintf(iw, " Ken Bagchi \\ / Version 0.0\n");
fprintf(iw, " \n");
}
@ -64,7 +64,7 @@ void qs_header(FILE *iw) {
fprintf(iw, " *** \n");
fprintf(iw, " \n");
fprintf(iw, " MPI Festkoerperforschung Stuttgart \n");
fprintf(iw, " 1999 \n");
fprintf(iw, " 2024 \n");
fprintf(iw, " \n");
fprintf(iw, " Version 0.0 \n");
fprintf(iw, " \n");
@ -87,7 +87,7 @@ void wave_header(FILE *iw) {
fprintf(iw, " **** **** **** **** ****** ******** \n");
fprintf(iw, " \n");
fprintf(iw, " MPI Festkoerperforschung Stuttgart \n");
fprintf(iw, " 1999 \n");
fprintf(iw, " 2024 \n");
fprintf(iw, " \n");
fprintf(iw, " Version 0.0 \n");
fprintf(iw, " \n");
@ -110,7 +110,7 @@ void faust_header(FILE *iw) {
fprintf(iw, " **** **** **** **** ****** ******** \n");
fprintf(iw, " \n");
fprintf(iw, " MPI Festkoerperforschung Stuttgart \n");
fprintf(iw, " 1999 \n");
fprintf(iw, " 2024 \n");
fprintf(iw, " \n");
fprintf(iw, " Version 0.0 \n");
fprintf(iw, " \n");

0
src/initialize_ewalds.F → src/initialize_ewalds.c
View file

0
src/initialize_extended_types.F → src/initialize_extended_types.c
View file

0
src/initialize_molecule_types.F → src/initialize_molecule_types.c
View file

0
src/initialize_particle_types.F → src/initialize_particle_types.c
View file

0
src/input_types.F → src/input_types.c
View file

0
src/integrator.F → src/integrator.c
View file

0
src/io_parameters.F → src/io_parameters.c
View file

2
src/kinds.cpp
View file

@ -1,7 +1,7 @@
/*
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
! Copyright (C) 2024 CP2K developers group !
!-----------------------------------------------------------------------------!
Data type definitions; tested on:
- IBM AIX xlf90

0
src/lapack.F → src/lapack.c
View file

32
src/linklist.h Normal file
View file

@ -0,0 +1,32 @@
#ifndef _LINKLIST_H
#define _LINKLIST_H
#ifdef __cplusplus
namespace CP2K_NS {
#endif
typedef struct cell_atoms_linklist_type {
int atom;
bool used;
cell_atoms_linklist_type *next;
} cell_atoms_linklist_type;
typedef struct cell_atoms_type {
cell_atoms_linklist_type *ll;
} cell_atoms_type;
typedef struct neighbor_cells_type {
int **index;
} neighbor_cells_type;
typedef struct cell_neighborlist_type {
int *natoms;
neighbor_cells_type **neighbor;
cell_atoms_type *ptype;
} cell_neighborlist_type;
#ifdef __cplusplus
}
#endif
#endif

0
src/linklist_cell_list.F → src/linklist_cell_list.c
View file

0
src/linklist_cell_types.F → src/linklist_cell_types.c
View file

0
src/linklist_control.F → src/linklist_control.c
View file

0
src/linklist_utilities.F → src/linklist_utilities.c
View file

0
src/linklist_verlet_list.F → src/linklist_verlet_list.c
View file

0
src/linklists.F → src/linklists.c
View file

0
src/machine.F → src/machine.c
View file

0
src/machine_absoft.F → src/machine_absoft.c
View file

0
src/machine_aix.F → src/machine_aix.c
View file

0
src/machine_dec.F → src/machine_dec.c
View file

0
src/machine_irix.F → src/machine_irix.c
View file

0
src/machine_pgi.F → src/machine_pgi.c
View file

0
src/machine_t3e.F → src/machine_t3e.c
View file

0
src/mathconstants.F → src/mathconstants.c
View file

499
src/md.F
View file

@ -1,499 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE md
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, stop_parser
USE particle_types, ONLY : particle_type
USE string_utilities, ONLY : uppercase, xstring
USE util, ONLY : gasdev
IMPLICIT NONE
PRIVATE
PUBLIC :: read_md_section, simulation_parameters_type, &
initialize_velocities, thermodynamic_type, mdio_parameters_type
TYPE simulation_parameters_type
INTEGER :: nsteps
REAL ( dbl ) :: dt
REAL ( dbl ) :: temp_ext
REAL ( dbl ) :: p_ext
CHARACTER ( LEN = 20 ) :: ensemble
CHARACTER ( LEN = 20 ) :: read_type
CHARACTER ( LEN = 20 ) :: program
LOGICAL :: constraint
INTEGER :: nfree
INTEGER :: nc
INTEGER :: nyosh
INTEGER :: nhclen
INTEGER, DIMENSION ( 3 ) :: n_cell
REAL ( dbl ) :: tau_nhc
REAL ( dbl ) :: tau_cell
REAL ( dbl ), POINTER, DIMENSION ( : ) :: dt_yosh
REAL ( dbl ) :: shake_tol
REAL ( dbl ) :: verlet_skin
END TYPE simulation_parameters_type
TYPE thermodynamic_type
REAL ( dbl ), DIMENSION ( 3, 3 ) :: ptens, pv_kin, pv, pv_const
REAL ( dbl ) :: kin, pot, gspace, e_self, e_neut, e_bonded
END TYPE thermodynamic_type
TYPE mdio_parameters_type
INTEGER :: icrd, ivel, iptens, iener, itemp, idump, iscreen
CHARACTER ( LEN = 50 ) :: crd_file_name
CHARACTER ( LEN = 50 ) :: vel_file_name
CHARACTER ( LEN = 50 ) :: ptens_file_name
CHARACTER ( LEN = 50 ) :: ener_file_name
CHARACTER ( LEN = 50 ) :: temp_file_name
CHARACTER ( LEN = 50 ) :: dump_file_name
END TYPE mdio_parameters_type
CONTAINS
!!>---------------------------------------------------------------------------!
!! SECTION: &md ... &end !
!! !
!! ensemble: [nve,nvt] ensemble type !
!! steps: n total number of steps !
!! timestep: dt time step [fs] !
!! temperature: T temperature [K] !
!! restart: [pos,posvel] start type !
!! constraints: [on,off] constraints !
!! shake: eps convergence for Shake algorithm !
!! verlet_skin vs buffer size for Verlet lists !
!! cells x y z Number of cells in each direction !
!! Nose_parameter & parameter of Nose-Hoover-chains !
!! length: nhclen & length of thermostat !
!! Yoshida: nyosh & order of Yoshida integrator !
!! mts: nc & number of multiple time steps !
!! timecon: tau_nhc NHC time constant !
!! Barostat_parameter parameter of barostat !
!! pressure: P & pressure [bar] !
!! timecon: tau_cell & barostat time constant !
!! FILES & !
!! filebody "filename" & !
!! coordinates "filename.crd" & !
!! velocities "filename.vel" & !
!! pressure "filename.ptens" & !
!! energies "filename.ener" & !
!! temperature "filename.temp" & !
!! dump "filename.dum" !
!! PRINT & !
!! screen iscreen & !
!! files ifiles & !
!! coordinates icrd & !
!! velocities ivel & !
!! pressure iptens & !
!! energies iener & !
!! temperature itemp & !
!! dump idump !
!! !
!! !
!!<---------------------------------------------------------------------------!
SUBROUTINE read_md_section ( simpar, mdpar, mdio )
IMPLICIT NONE
! Arguments
TYPE ( simulation_parameters_type ), INTENT ( OUT ) :: simpar
TYPE ( global_environment_type ), INTENT ( IN ) :: mdpar
TYPE ( mdio_parameters_type ), INTENT ( OUT ) :: mdio
! Locals
INTEGER :: ierror, ilen, msglen, ia, ie, iw, allgrp, source
CHARACTER ( LEN = 20 ) :: string
CHARACTER ( LEN = 5 ) :: label
!------------------------------------------------------------------------------
!..defaults
simpar % ensemble = 'NVE'
simpar % nsteps = 0
simpar % dt = 0.0_dbl
simpar % temp_ext = 0.0_dbl
simpar % p_ext = 0.0_dbl
simpar % read_type = 'POS'
simpar % constraint = .FALSE.
simpar % shake_tol = 1.0E-6_dbl
simpar % nhclen = 1
simpar % nc = 1
simpar % n_cell = 1
simpar % nyosh = 1
simpar % tau_nhc = 1000.0_dbl
simpar % tau_cell = 1000.0_dbl
simpar % verlet_skin = 1.0_dbl
iw = mdpar % scr
!..filenames
CALL xstring ( mdpar % project_name, ia, ie )
mdio % crd_file_name = mdpar % project_name(ia:ie) // '.crd'
mdio % vel_file_name = mdpar % project_name(ia:ie) // '.vel'
mdio % ptens_file_name = mdpar % project_name(ia:ie) // '.ptens'
mdio % ener_file_name = mdpar % project_name(ia:ie) // '.ener'
mdio % temp_file_name = mdpar % project_name(ia:ie) // '.temp'
mdio % dump_file_name = mdpar % project_name(ia:ie) // '.dump'
!..print frequency
mdio % icrd = 1
mdio % ivel = 1
mdio % iptens = 1
mdio % iener = 1
mdio % itemp = 1
mdio % idump = 1
mdio % iscreen = 1
!..parse the input section
label = '&MD'
CALL parser_init(mdpar % input_file_name,label,ierror,mdpar)
IF (ierror /= 0 ) THEN
IF( mdpar % ionode ) &
WRITE ( iw, '( a )' ) ' No input section &MD found '
ELSE
CALL read_line
DO WHILE ( test_next() /= 'X' )
ilen = 6
CALL cfield ( string, ilen )
CALL uppercase ( string )
SELECT CASE ( string )
CASE DEFAULT
CALL p_error()
CALL stop_parser ( 'read_md_section', 'unknown option' )
CASE ( 'ENSEMB' )
ilen = 0
CALL cfield ( string, ilen )
CALL uppercase ( string )
simpar % ensemble = string ( 1:ilen )
CASE ( 'STEPS' )
simpar % nsteps = get_int()
CASE ( 'VERLET' )
simpar % verlet_skin = get_real()
CASE ( 'CELL' )
simpar % n_cell(1) = get_int()
simpar % n_cell(2) = get_int()
simpar % n_cell(3) = get_int()
CASE ( 'TIMEST' )
simpar % dt = get_real()
CASE ( 'TEMPER' )
simpar % temp_ext = get_real()
CASE ( 'RESTAR' )
ilen = 0
CALL cfield ( string, ilen )
CALL uppercase ( string )
simpar % read_type = string ( 1:ilen )
CASE ( 'CONSTR' )
ilen = 0
CALL cfield ( string, ilen )
CALL uppercase ( string )
IF ( string(1:2)=='ON' ) simpar % constraint = .TRUE.
IF ( string(1:3)=='OFF' ) simpar % constraint = .FALSE.
CASE ( 'SHAKE' )
simpar % shake_tol = get_real()
CASE ( 'NOSE_P' )
DO WHILE ( test_next() == 'C' )
ilen = 6
CALL cfield ( string, ilen )
CALL uppercase ( string )
SELECT CASE ( string )
CASE DEFAULT
CALL p_error()
CALL stop_parser ( 'read_md_section%nose', &
'unknown suboption' )
CASE ( 'LENGTH' )
simpar % nhclen = get_int()
CASE ( 'YOSHID' )
simpar % nyosh = get_int()
CASE ( 'TIMECO' )
simpar % tau_nhc = get_real()
CASE ( 'MTS' )
simpar % nc = get_int()
END SELECT
END DO
CASE ( 'BAROST' )
DO WHILE ( test_next() == 'C' )
ilen = 6
CALL cfield ( string, ilen )
CALL uppercase ( string )
SELECT CASE ( string )
CASE DEFAULT
CALL p_error()
CALL stop_parser ( 'read_md_section%barost', &
'unknown suboption' )
CASE ( 'TIMECO' )
simpar % tau_cell = get_real()
CASE ( 'PRESSU' )
simpar % p_ext = get_real()
END SELECT
END DO
CASE ( 'FILES' )
DO WHILE ( test_next() == 'C' )
ilen = 6
CALL cfield ( string, ilen )
CALL uppercase ( string )
SELECT CASE ( string )
CASE DEFAULT
CALL p_error()
CALL stop_parser ( 'read_md_section%files', &
'unknown suboption' )
CASE ( 'FILEBO' )
ilen = 0
CALL cfield ( string, ilen )
mdio % crd_file_name = string(1:ilen) // '.crd'
mdio % vel_file_name = string(1:ilen) // '.vel'
mdio % ptens_file_name = string(1:ilen) // '.ptens'
mdio % ener_file_name = string(1:ilen) // '.ener'
mdio % temp_file_name = string(1:ilen) // '.temp'
mdio % dump_file_name = string(1:ilen) // '.dump'
CASE ( 'COORDI' )
ilen = 0
CALL cfield ( mdio % crd_file_name, ilen )
CASE ( 'ENERGI' )
ilen = 0
CALL cfield ( mdio % ener_file_name, ilen )
CASE ( 'PRESSU' )
ilen = 0
CALL cfield ( mdio % ptens_file_name, ilen )
CASE ( 'VELOCI' )
ilen = 0
CALL cfield ( mdio % vel_file_name, ilen )
CASE ( 'TEMPER' )
ilen = 0
CALL cfield ( mdio % temp_file_name, ilen )
CASE ( 'DUMP' )
ilen = 0
CALL cfield ( mdio % dump_file_name, ilen )
END SELECT
END DO
CASE ( 'PRINT' )
DO WHILE ( test_next() == 'C' )
ilen = 6
CALL cfield ( string, ilen )
CALL uppercase ( string )
SELECT CASE ( string )
CASE DEFAULT
CALL p_error()
CALL stop_parser ( 'read_md_section%print', &
'unknown suboption' )
CASE ( 'SCREEN' )
mdio % iscreen = get_int()
CASE ( 'FILES' )
ia = get_int()
mdio % icrd = ia
mdio % ivel = ia
mdio % iener = ia
mdio % iptens = ia
mdio % itemp = ia
mdio % idump = ia
CASE ( 'COORDI' )
mdio % icrd = get_int()
CASE ( 'ENERGI' )
mdio % iener = get_int()
CASE ( 'PRESSU' )
mdio % iptens = get_int()
CASE ( 'VELOCI' )
mdio % ivel = get_int()
CASE ( 'TEMPER' )
mdio % itemp = get_int()
CASE ( 'DUMP' )
mdio % idump = get_int()
END SELECT
END DO
END SELECT
CALL read_line
END DO
END IF
CALL parser_end
!..end of parsing the input section
!..write some information to output
IF (mdpar % ionode .AND. mdpar % print_level>=0) THEN
WRITE ( iw, '( A )' ) ' MD| Molecular Dynamics Protocol '
WRITE ( iw, '( A,T61,A )' ) ' MD| Ensemble Type ', &
ADJUSTR ( simpar % ensemble )
WRITE ( iw, '( A,T61,A )' ) ' MD| Restart Type ', &
ADJUSTR ( simpar % read_type )
WRITE ( iw, '( A,T71,I10 )' ) ' MD| Number of Time Steps ', &
simpar % nsteps
WRITE ( iw, '( A,T71,F10.2 )' ) ' MD| Time Step [fs] ', simpar % dt
WRITE ( iw, '( A,T71,F10.2 )' ) ' MD| Temperature [K] ', &
simpar % temp_ext
WRITE ( iw, '( A,T71,F10.2 )' ) &
' MD| Verlet buffer skin [Angstrom] ', simpar % verlet_skin
IF ( ( simpar % ensemble ( 1:3 ) == 'NVT' ) &
.OR. (simpar % ensemble ( 1:3 ) == 'NPT' ) ) THEN
WRITE ( iw, '( A )' ) ' MD| Nose-Hoover-Chain parameters '
WRITE ( iw, '( A,T77,I4 )' ) ' MD| Nose-Hoover-Chain length ', &
simpar % nhclen
WRITE ( iw, '( A,T71,F10.4 )' ) &
' MD| Nose-Hoover-Chain time constant [fs] ', simpar % tau_nhc
WRITE ( iw, '( A,T77,I4 )' ) ' MD| Order of Yoshida integrator ', &
simpar % nyosh
WRITE ( iw, '( A,T77,I4 )' ) ' MD| Number of multiple time steps', &
simpar % nc
END IF
IF ( simpar % ensemble ( 1:3 ) == 'NPT' ) THEN
WRITE ( iw, '( A, T71, F10.2 )' ) &
' MD| Pressure [bar] ', simpar % p_ext
WRITE ( iw, '( A, T71, F10.4 )' ) &
' MD| Barostat time constant [fs] ', simpar % tau_cell
END IF
IF ( simpar % constraint ) THEN
WRITE ( iw, '( A )' ) ' MD| Constraints activated '
WRITE ( iw, '( A,T71,G10.4 )' ) ' MD| Tolerance for shake ', &
simpar % shake_tol
END IF
WRITE ( iw, '( A,T63,i10,A )' ) ' MD| Print MD information every', &
mdio % iscreen, ' step(s)'
WRITE ( iw, '( A,T20,A,T71,A )' ) ' MD| File type', &
'Print frequency[steps]', 'File names'
WRITE ( iw, '( A,T25,i5,T31,A )' ) ' MD| Coordinates', mdio % icrd, &
ADJUSTR ( mdio % crd_file_name ( 1:50 ) )
WRITE ( iw, '( A,T25,i5,T31,A )' ) ' MD| Velocities', mdio % ivel, &
ADJUSTR ( mdio % vel_file_name ( 1:50 ) )
WRITE ( iw, '( A,T25,i5,T31,A )' ) ' MD| Energies', mdio % iener, &
ADJUSTR ( mdio % ener_file_name ( 1:50 ) )
WRITE ( iw, '( A,T25,i5,T31,A )' ) ' MD| Temperature', mdio % itemp, &
ADJUSTR ( mdio % temp_file_name ( 1:50 ) )
WRITE ( iw, '( A,T25,i5,T31,A )' ) ' MD| Pressure Tensor', mdio % iptens, &
ADJUSTR ( mdio % ptens_file_name ( 1:50 ) )
WRITE ( iw, '( A,T25,i5,T31,A )' ) ' MD| Dump', mdio % idump, &
ADJUSTR ( mdio % dump_file_name ( 1:50 ) )
WRITE (iw,*)
END IF
END SUBROUTINE read_md_section
!******************************************************************************
! Initializes the velocities to the Maxwellian distribution
SUBROUTINE initialize_velocities ( simpar, part, mdpar )
IMPLICIT NONE
! Arguments
TYPE ( simulation_parameters_type ), INTENT ( IN ) :: simpar
TYPE ( particle_type ), INTENT ( INOUT ), DIMENSION ( : ) :: part
TYPE ( global_environment_type ), INTENT ( INOUT ) :: mdpar
! Locals
INTEGER :: i, id, natoms, iw
REAL ( dbl ), DIMENSION ( 3 ) :: vcom
REAL ( dbl ), DIMENSION ( 3 ) :: sum
REAL ( dbl ) :: denom, sc, mass
REAL ( dbl ) :: akin, temp
LOGICAL :: domdec
!------------------------------------------------------------------------------
iw = mdpar % scr
IF (mdpar % ionode .AND. mdpar % print_level>1) THEN
WRITE ( iw, '( A, A, A )' ) ' **********************', &
' begin of velocity initialization ', '***********************'
END IF
! Initializing parameters
sum = 0.0_dbl
vcom = 0.0_dbl
denom = 0.0_dbl
natoms = SIZE ( part )
! initializing velocities
IF (mdpar % ionode) THEN
DO i = 1, natoms
part(i) % v(1) = gasdev(mdpar % idum)
part(i) % v(2) = gasdev(mdpar % idum)
part(i) % v(3) = gasdev(mdpar % idum)
END DO
akin = 0.0_dbl
DO i = 1, natoms
mass = part(i) % prop % mass
akin = akin + 0.5_dbl*(mass*part(i) % v(1)*part(i) % v(1))
akin = akin + 0.5_dbl*(mass*part(i) % v(2)*part(i) % v(2))
akin = akin + 0.5_dbl*(mass*part(i) % v(3)*part(i) % v(3))
END DO
temp = 2.0_dbl * akin / REAL ( simpar % nfree, dbl )
! scale velocities to get the correct initial temperature
sc = SQRT ( simpar % temp_ext / temp )
DO i = 1, natoms
part(i) % v(1) = sc*part(i) % v(1)
part(i) % v(2) = sc*part(i) % v(2)
part(i) % v(3) = sc*part(i) % v(3)
END DO
! Computing the COM velocities
DO i = 1, natoms
mass = part(i) % prop % mass
vcom(1) = vcom(1) + mass*part(i) % v(1)
vcom(2) = vcom(2) + mass*part(i) % v(2)
vcom(3) = vcom(3) + mass*part(i) % v(3)
denom = denom + mass
END DO
! Moving the velocities so that the COM velocity is zero
vcom = vcom / denom
DO i = 1, natoms
mass = part(i) % prop % mass
part(i) % v ( : ) = part(i) % v ( : ) - vcom ( : )
sum ( : ) = sum ( : ) + mass*part(i) % v ( : )
END DO
END IF
#if defined(__parallel)
DO i = 1, natoms
CALL mp_bcast(part(i) % v,mdpar % source,mdpar % group)
END DO
#endif
IF ( mdpar % ionode .AND. mdpar % print_level > 1 ) THEN
WRITE ( iw, '( A, T69, E12.6 )' ) &
' Centre of mass velocity in direction x:', sum ( 1 )
WRITE ( iw, '( A, T69, E12.6 )' ) &
' Centre of mass velocity in direction y:', sum ( 2 )
WRITE ( iw, '( A, T69, E12.6 )' ) &
' Centre of mass velocity in direction z:', sum ( 3 )
WRITE ( iw, '( A, A, A, / )' ) ' ***********************', &
' end of velocity initialization ', '************************'
END IF
END SUBROUTINE initialize_velocities
!******************************************************************************
END MODULE md

6
src/md.h Normal file
View file

@ -0,0 +1,6 @@
#ifndef _MD_H
#define _MD_H
#endif

0
src/mol_force.F → src/mol_force.c
View file

165
src/molecule.h Normal file
View file

@ -0,0 +1,165 @@
#ifndef _MOLECULE_H
#define _MOLECULE_H
#include "particle.h"
#include "linklist.h"
#ifdef __cplusplus
namespace CP2K_NS {
#endif
typedef struct linklist_images {
double vec[3];
linklist_images *next;
} linklist_images;
typedef struct linklist_exclusion {
particle_type *p;
linklist_exclusion *next;
} linklist_exclusion;
typedef struct localization_list {
bool all;
int nloc;
int *atomlist;
} localization_list;
typedef struct linklist_neighbor {
particle_type *p;
linklist_images *image;
int index, nimages;
linklist_neighbor *next;
} linklist_neighbor;
typedef struct particle_node_type {
particle_type *p;
int in_mol, iat_in_mol;
int nneighbor, nexcl, nsneighbor;
linklist_neighbor *nl;
linklist_neighbor *sl;
linklist_exclusion *ex;
localization_list *llist;
} particle_node_type;
typedef struct linklist_atoms {
particle_node_type *part;
linklist_atoms *next;
} linklist_atoms;
typedef struct bond_parameters_type {
char type[21];
double r0, k;
} bond_parameters_type;
typedef struct linklist_bonds {
particle_type *p1, *p2;
bond_parameters_type *bond_param;
linklist_bonds *next;
int index[2];
} linklist_bonds;
typedef struct bend_parameters_type {
char type[21];
double theta0, k;
} bend_parameters_type;
typedef struct linklist_bends {
particle_type *p1, *p2, *p3;
bend_parameters_type *bend_param;
linklist_bends *next;
int index[3];
} linklist_bends;
typedef struct torsion_parameters_type {
char type[21];
double parm[7];
} torsion_parameters_type;
typedef struct linklist_torsions {
particle_type *p1, *p2, *p3, *p4;
torsion_parameters_type *torsion_param;
linklist_torsions *next;
int index[4];
} linklist_torsions;
typedef struct linklist_dist_constraints {
particle_type *p1, *p2;
double *distance;
double del_lam, lambda, sigma;
double rold1[3], rold2[3], fc[3];
double scale, scale_old;
linklist_dist_constraints *next;
} linklist_dist_constraints;
typedef struct linklist_g3x3_constraints {
particle_type *p1, *p2, *p3;
double *d12, *d13, *d23;
double fc1[3], fc2[3], fc3[3];
double lambda[3], del_lam[3];
double rold1[3], rold2[3], rold3[3];
double vold1[3], vold2[3], vold3[3];
double v1[3], v2[3], v3[3];
double scale, scale_old;
linklist_g3x3_constraints *next;
} linklist_g3x3_constraints;
typedef struct intra_parameters_type {
bond_parameters_type **bond_param;
bend_parameters_type ***bend_param;
torsion_parameters_type ***torsion_param;
double **constraint_distance;
} intra_parameters_type;
typedef struct molecule_parameter_type {
char name[21];
int natom;
int natom_type;
char *aname[21];
char *alabel[21];
double *aweight;
double *acharge;
int nbonds;
int **bonds;
int nbends;
int **bends;
int ntorsions ;
int **torsions;
int ndcon ;
int **ndc;
int n3x3con;
int **n3x3c;
} molecule_parameter_type;
typedef struct molecule_type {
char file_name[21];
int num_mol;
char nhcopt[21];
molecule_parameter_type molpar;
} molecule_type;
typedef struct molecule_structure_type {
char mol_name[31]; // name of molecule
int moltype; // pointer to molecular type
int imol; // position counter
int atombase; // position counter of first atom in molecule
int natoms_mol; // # of atoms in molecule
int nvt_num; // thermostat number that it is connected.
// Used for constraints
linklist_atoms *ll_atoms;
int nbonds_mol; // # of bonds in molecule
linklist_bonds *ll_bonds;
int nbends_mol; // # of bends in molecule
linklist_bends *ll_bends;
int ntorsions_mol; // # of torsions in molecule
linklist_torsions *ll_torsions;
int ndistconst_mol; // # of distance constraints in molecule
linklist_dist_constraints *ll_dist_const;
int ng3x3const_mol; // # of 3x3 group constraints in molecule
linklist_g3x3_constraints *ll_g3x3_const;
} molecule_structure_type;
#ifdef __cplusplus
}
#endif
#endif

0
src/molecule_input.F → src/molecule_input.c
View file

0
src/molecule_types.F → src/molecule_types.c
View file

339
src/mp.F → src/mp.cpp
View file

@ -1,218 +1,127 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
/*---------------------------------------------------------------------------*/
/* CP2K: A general program to perform molecular dynamics simulations */
/* Copyright (C) 2024 CP2K developers group */
/*---------------------------------------------------------------------------*/
MODULE mp
USE kinds, ONLY : dbl
USE parallel_include
PRIVATE
PUBLIC :: mp_start, mp_end, mp_environ, mp_group, mp_cart_create, &
mp_bcast, mp_stop, mp_sum, mp_max, mp_min, mp_sync
INTERFACE mp_bcast
MODULE PROCEDURE mp_bcast_i1, mp_bcast_r1, mp_bcast_c1, mp_bcast_z, &
mp_bcast_iv, mp_bcast_rv, mp_bcast_cv, mp_bcast_l, mp_bcast_rm, &
mp_bcast_cm, mp_bcast_im
END INTERFACE
INTERFACE mp_sum
MODULE PROCEDURE mp_sum_i1, mp_sum_r1, mp_sum_c1, mp_sum_iv, &
mp_sum_rv, mp_sum_cv, mp_sum_rm, mp_sum_cm
END INTERFACE
INTERFACE mp_max
MODULE PROCEDURE mp_max_i, mp_max_r, mp_max_rv
END INTERFACE
INTERFACE mp_min
MODULE PROCEDURE mp_min_i, mp_min_r, mp_min_rv
END INTERFACE
CONTAINS
#include <iostream>
!******************************************************************************
!..mp_start
SUBROUTINE mp_start()
IMPLICIT NONE
! Locals
INTEGER :: ierr
ierr = 0
#if defined(__parallel)
CALL mpi_init ( ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8000 )
#endif
END SUBROUTINE mp_start
void mp_start() {
!******************************************************************************
int ierr = 0;
#if defined(__parallel)
mpi_init(ierr);
if (ierr != 0) mp_stop(8000);
#endif
!..mp_end
SUBROUTINE mp_end()
IMPLICIT NONE
! Locals
INTEGER :: ierr
ierr = 0
#if defined(__parallel)
CALL mpi_finalize ( ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8999 )
#endif
END SUBROUTINE mp_end
}
!******************************************************************************
void mp_end() {
!..mp_stop
SUBROUTINE mp_stop ( code )
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( IN ) :: code
! Locals
INTEGER :: numtask, taskid, groupid
#if defined(__parallel)
CALL mpi_abort ( mpi_comm_world, code )
#endif
CALL mp_environ ( numtask, taskid, groupid )
WRITE(*,'(/,A,T71,I10,/)') ' CP2K| Stopped by processor number',taskid
int ierr = 0;
#if defined(__parallel)
mpi_finalize(ierr);
if (ierr != 0) mp_stop(8999);
#endif
STOP "mp_stop"
}
END SUBROUTINE mp_stop
void mp_stop(int code) {
!******************************************************************************
int numtask, taskid, groupid;
!..mp_sync
SUBROUTINE mp_sync ( group )
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( IN ) :: group
! Locals
INTEGER :: ierr
#if defined(__parallel)
CALL mpi_barrier ( group, ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8400 )
#endif
END SUBROUTINE mp_sync
#if defined(__parallel)
mpi_abort(mpi_comm_world, code);
#endif
!******************************************************************************
mp_environ(numtask, taskid, groupid);
!..mp_environ
SUBROUTINE mp_environ ( numtask, taskid, groupid )
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( OUT ) :: numtask, taskid, groupid
! Locals
INTEGER :: ierr
ierr = 0
numtask = 1
taskid = 0
groupid = 0
#if defined(__parallel)
CALL mpi_comm_rank ( mpi_comm_world, taskid, ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8001 )
CALL mpi_comm_size ( mpi_comm_world, numtask, ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8002 )
groupid = mpi_comm_world
#endif
END SUBROUTINE mp_environ
std::cout << " CP2K| Stopped by processor number " << taskid;
!******************************************************************************
std::cout << "mp_stop" << std::endl;
exit(1);
}
!..mp_group
SUBROUTINE mp_group(group_list,group_size,base_group,groupid)
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( IN ) :: group_list ( : ), group_size, base_group
INTEGER, INTENT ( OUT ) :: groupid
! Locals
INTEGER :: base, newgroup, ierr
ierr = 0
groupid = base_group
#if defined(__parallel)
CALL mpi_comm_group ( base_group, base, ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8010 )
CALL mpi_group_incl ( base, group_size, group_list, newgroup, ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8011 )
CALL mpi_comm_create ( base_group, newgroup, groupid, ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8012 )
#endif
END SUBROUTINE mp_group
void mp_sync(int group) {
!******************************************************************************
int ierr = 0;
#if defined(__parallel)
mpi_barrier(group, ierr);
if (ierr != 0) mp_stop(8400);
#endif
}
!..mp_cart_create
SUBROUTINE mp_cart_create ( comm_old, ndims, dims, pos, comm_cart )
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( IN ) :: comm_old, ndims
INTEGER, INTENT ( OUT ) :: dims ( : ), pos ( : ), comm_cart
! Locals
INTEGER :: ierr, nodes
LOGICAL :: period(1:ndims), reorder
ierr = 0
dims ( 1:ndims ) = 1
pos ( 1:ndims ) = 1
comm_cart = comm_old
#if defined(__parallel)
dims ( 1:ndims ) = 0
CALL mpi_comm_size ( comm_old, nodes, ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8020 )
CALL mpi_dims_create ( nodes, ndims, dims, ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8021 )
reorder = .TRUE.
period = .TRUE.
CALL mpi_cart_create ( comm_old, ndims, dims, period, reorder, comm_cart, &
ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8022 )
CALL mpi_cart_get ( comm_cart, ndims, dims, period, pos, ierr )
IF ( ierr /= 0 ) CALL mp_stop ( 8023 )
#endif
END SUBROUTINE mp_cart_create
void mp_environ(int numtask, int taskid, int groupid) {
!******************************************************************************
int ierr = 0;
numtask = 1;
taskid = 0;
groupid = 0;
#if defined(__parallel)
mpi_comm_rank(mpi_comm_world, taskid, ierr);
if (ierr != 0) mp_stop(8001);
mpi_comm_size(mpi_comm_world, numtask, ierr);
if (ierr != 0) mp_stop(8002);
groupid = mpi_comm_world;
#endif
}
void mp_group(int *group_list, int group_size, int base_group, int *groupid) {
int base, newgroup, ierr;
ierr = 0;
*groupid = base_group;
#if defined(__parallel)
mpi_comm_group(base_group, base, ierr);
if (ierr != 0) mp_stop(8010);
mpi_group_incl(base, group_size, group_list, newgroup, ierr);
if (ierr != 0) mp_stop(8011);
mpi_comm_create(base_group, newgroup, groupid, ierr);
if (ierr != 0) mp_stop(8012);
#endif
}
void mp_cart_create (int comm_old, int ndims, int *dims, int *pos, int *comm_cart ) {
size_t i;
int ierr, nodes;
bool period[ndims], reorder;
ierr = 0;
for (i = 0; i < ndims; i++) {
dims[i] = 1;
pos[i] = 1;
}
comm_cart = comm_old;
#if defined(__parallel)
for (i = 0; i < ndims; i++) {
dims[i] = 0;
}
mpi_comm_siz ( comm_old, nodes, ierr );
if ( ierr != 0 ) mp_stop ( 8020 );
mpi_dims_create ( nodes, ndims, dims, ierr );
if ( ierr != 0 ) mp_stop ( 8021 );
reorder = true;
period = true;
mpi_cart_create ( comm_old, ndims, dims, period, reorder, comm_cart, ierr );
if ( ierr != 0 ) mp_stop ( 8022 );
mpi_cart_get ( comm_cart, ndims, dims, period, pos, ierr );
IF ( ierr != 0 ) mp_stop ( 8023 );
#endif
}
!..mp_bcast
SUBROUTINE mp_bcast_i1(msg,source,gid)
@ -220,7 +129,7 @@ SUBROUTINE mp_bcast_i1(msg,source,gid)
INTEGER :: msg
INTEGER :: source, gid
INTEGER :: msglen, ierr
ierr = 0
msglen = 1
#if defined(__parallel)
@ -233,7 +142,7 @@ SUBROUTINE mp_bcast_iv(msg,source,gid)
INTEGER :: msg ( : )
INTEGER :: source, gid
INTEGER :: msglen, ierr
ierr = 0
msglen = size(msg)
#if defined(__parallel)
@ -246,7 +155,7 @@ SUBROUTINE mp_bcast_im(msg,source,gid)
INTEGER :: msg ( :, : )
INTEGER :: source, gid
INTEGER :: msglen, ierr
ierr = 0
msglen = size(msg)
#if defined(__parallel)
@ -259,7 +168,7 @@ SUBROUTINE mp_bcast_r1(msg,source,gid)
REAL ( dbl ) :: msg
INTEGER :: msglen, source, gid
INTEGER :: ierr
ierr = 0
msglen = 1
#if defined(__parallel)
@ -272,7 +181,7 @@ SUBROUTINE mp_bcast_rv(msg,source,gid)
REAL ( dbl ) :: msg ( : )
INTEGER :: source, gid
INTEGER :: msglen, ierr
ierr = 0
msglen = size(msg)
#if defined(__parallel)
@ -285,7 +194,7 @@ SUBROUTINE mp_bcast_rm(msg,source,gid)
REAL ( dbl ) :: msg ( :, : )
INTEGER :: source, gid
INTEGER :: msglen, ierr
ierr = 0
msglen = size(msg)
#if defined(__parallel)
@ -298,7 +207,7 @@ SUBROUTINE mp_bcast_c1(msg,source,gid)
COMPLEX ( dbl ) :: msg
INTEGER :: source, gid
INTEGER :: msglen, ierr
ierr = 0
msglen = 1
#if defined(__parallel)
@ -311,7 +220,7 @@ SUBROUTINE mp_bcast_cv(msg,source,gid)
COMPLEX ( dbl ) :: msg ( : )
INTEGER :: source, gid
INTEGER :: msglen, ierr
ierr = 0
msglen = size(msg)
#if defined(__parallel)
@ -324,7 +233,7 @@ SUBROUTINE mp_bcast_cm(msg,source,gid)
COMPLEX ( dbl ) :: msg ( :, : )
INTEGER :: source, gid
INTEGER :: msglen, ierr
ierr = 0
msglen = size(msg)
#if defined(__parallel)
@ -337,7 +246,7 @@ SUBROUTINE mp_bcast_l(msg,source,gid)
LOGICAL :: msg
INTEGER :: source, gid
INTEGER :: msglen, ierr
ierr = 0
msglen = 1
#if defined(__parallel)
@ -351,7 +260,7 @@ SUBROUTINE mp_bcast_z(msg,source,gid)
INTEGER :: source, gid
INTEGER :: msglen, ierr, i
INTEGER, ALLOCATABLE :: imsg ( : )
ierr = 0
msglen = len(msg)
#if defined(__parallel)
@ -381,7 +290,7 @@ SUBROUTINE mp_sum_i1(msg,gid)
INTEGER, INTENT ( INOUT ) :: msg
INTEGER, INTENT ( IN ) :: gid
INTEGER :: msglen, res, ierr
msglen = 1
#if defined(__parallel)
CALL mpi_allreduce(msg,res,msglen,mpi_integer,mpi_sum,gid,ierr)
@ -411,7 +320,7 @@ SUBROUTINE mp_sum_r1(msg,gid)
INTEGER, INTENT ( IN ) :: gid
INTEGER :: msglen, ierr
REAL ( dbl ) :: res
msglen = 1
#if defined(__parallel)
CALL mpi_allreduce(msg,res,msglen,mpi_double_precision,mpi_sum,gid, &
@ -462,7 +371,7 @@ SUBROUTINE mp_sum_c1(msg,gid)
INTEGER, INTENT ( IN ) :: gid
INTEGER :: msglen, ierr
COMPLEX ( dbl ) :: res
msglen = 2
#if defined(__parallel)
CALL mpi_allreduce(msg,res,msglen,mpi_double_precision,mpi_sum,gid, &
@ -600,8 +509,4 @@ SUBROUTINE mp_min_rv(msg,gid)
DEALLOCATE (res)
IF ( ierr /= 0 ) CALL mp_stop(8313)
#endif
END SUBROUTINE mp_min_rv
!******************************************************************************
END MODULE mp
END SUBROUTINE mp_min_rv

0
src/nhc.F → src/nhc.c
View file

0
src/nose.F → src/nose.c
View file

67
src/nrutil.cpp
View file

@ -1,10 +1,14 @@
#include <cmath>
#include <string>
#include <vector>
#include <algorithm>
namespace nrutil {
namespace CP2K_NS {
const double twopi = 6.283185307179586476925286766559005768394;
const int npar_arth = 16, npar2_arth = 8;
const int npar_geop = 4, npar2_geop = 2;
const int npar_arth = 16;
const int npar2_arth = 8;
const int npar_geop = 4;
const int npar2_geop = 2;
const int npar_cumsum = 16;
const int npar_cumprod = 8;
const int npar_poly = 8;
@ -232,7 +236,7 @@ namespace nrutil {
iminloc = imin(1)
END FUNCTION iminloc
SUBROUTINE assert1(n1,string)
void assert1(bool n1, std::string string) {
CHARACTER (len=*), INTENT (IN) :: string
LOGICAL, INTENT (IN) :: n1
@ -240,9 +244,9 @@ namespace nrutil {
WRITE (*,*) 'nrerror: an assertion failed with this tag:', string
STOP 'program terminated by assert1'
END IF
END SUBROUTINE assert1
}
SUBROUTINE assert2(n1,n2,string)
void assert2(bool n1, bool n2, std::string string) {
CHARACTER (len=*), INTENT (IN) :: string
LOGICAL, INTENT (IN) :: n1, n2
@ -250,9 +254,9 @@ namespace nrutil {
WRITE (*,*) 'nrerror: an assertion failed with this tag:', string
STOP 'program terminated by assert2'
END IF
END SUBROUTINE assert2
}
SUBROUTINE assert3(n1,n2,n3,string)
void assert3(bool n1, bool n2, bool n3, std::string string) {
CHARACTER (len=*), INTENT (IN) :: string
LOGICAL, INTENT (IN) :: n1, n2, n3
@ -260,9 +264,9 @@ namespace nrutil {
WRITE (*,*) 'nrerror: an assertion failed with this tag:', string
STOP 'program terminated by assert3'
END IF
END SUBROUTINE assert3
}
SUBROUTINE assert4(n1,n2,n3,n4,string)
void assert4(bool n1, bool n2, bool n3, bool n4, std::string string) {
CHARACTER (len=*), INTENT (IN) :: string
LOGICAL, INTENT (IN) :: n1, n2, n3, n4
@ -270,9 +274,9 @@ namespace nrutil {
WRITE (*,*) 'nrerror: an assertion failed with this tag:', string
STOP 'program terminated by assert4'
END IF
END SUBROUTINE assert4
}
SUBROUTINE assert_v(n,string)
void assert_v(bool n[], std::string string) {
CHARACTER (len=*), INTENT (IN) :: string
LOGICAL, DIMENSION (:), INTENT (IN) :: n
@ -280,7 +284,7 @@ namespace nrutil {
WRITE (*,*) 'nrerror: an assertion failed with this tag:', string
STOP 'program terminated by assert_v'
END IF
END SUBROUTINE assert_v
}
int assert_eq2(int n1, int n2, std::string string) {
@ -292,7 +296,7 @@ namespace nrutil {
}
}
FUNCTION assert_eq3(n1,n2,n3,string)
int assert_eq3(int n1, int n2, int n3, std::string) {
CHARACTER (len=*), INTENT (IN) :: string
INTEGER, INTENT (IN) :: n1, n2, n3
INTEGER :: assert_eq3
@ -303,9 +307,9 @@ namespace nrutil {
WRITE (*,*) 'nrerror: an assert_eq failed with this tag:', string
STOP 'program terminated by assert_eq3'
END IF
END FUNCTION assert_eq3
}
FUNCTION assert_eq4(n1,n2,n3,n4,string)
int assert_eq4(int n1, int n2, int n3, int n4, std::string) {
CHARACTER (len=*), INTENT (IN) :: string
INTEGER, INTENT (IN) :: n1, n2, n3, n4
INTEGER :: assert_eq4
@ -316,9 +320,9 @@ namespace nrutil {
WRITE (*,*) 'nrerror: an assert_eq failed with this tag:', string
STOP 'program terminated by assert_eq4'
END IF
END FUNCTION assert_eq4
}
FUNCTION assert_eqn(nn,string)
int assert_eqn(int nn[], std::string string) {
CHARACTER (len=*), INTENT (IN) :: string
INTEGER, DIMENSION (:), INTENT (IN) :: nn
INTEGER :: assert_eqn
@ -329,23 +333,22 @@ namespace nrutil {
WRITE (*,*) 'nrerror: an assert_eq failed with this tag:', string
STOP 'program terminated by assert_eqn'
END IF
END FUNCTION assert_eqn
}
SUBROUTINE nrerror(string)
CHARACTER (len=*), INTENT (IN) :: string
void nrerror(std::string string) {
printf("nrerror: %s", string);
printf("program terminated by nerror");
exit(1);
}
WRITE (*,*) 'nrerror: ', string
STOP 'program terminated by nrerror'
END SUBROUTINE nrerror
FUNCTION arth_r(first,increment,n)
REAL (sp), INTENT (IN) :: first, increment
INTEGER (i4b), INTENT (IN) :: n
REAL (sp), DIMENSION (n) :: arth_r
INTEGER (i4b) :: k, k2
REAL (sp) :: temp
double *arth_r(double first, double increment, int n) {
int k, k2;
double temp;
double *arth_r = new double[n];
if (n > 0) arth_r[0] = first;
IF (n>0) arth_r(1) = first
IF (n<=npar_arth) THEN
DO k = 2, n
arth_r(k) = arth_r(k-1) + increment
@ -364,7 +367,7 @@ namespace nrutil {
k = k2
END DO
END IF
END FUNCTION arth_r
}
FUNCTION arth_i(first,increment,n)
INTEGER (i4b), INTENT (IN) :: first, increment, n

0
src/output_utilities.F → src/output_utilities.c
View file

0
src/pair_potential.F → src/pair_potential.c
View file

0
src/parallel.F → src/parallel.c
View file

0
src/parallel_include.F → src/parallel_include.c
View file

168
src/parser.F → src/parser.cpp
View file

@ -1,117 +1,67 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
!-----------------------------------------------------------------------------!
!!> F R E E - F O R M A T I N P U T R O U T I N E S !
!! !
!! Function INSCAN(iunit,label) : looks for an input dection started by !
!! label un file associated with unit iunit !
!! Subroutine READ_LINE(iunit) : reads an input line, physical end of lines !
!! can be masked by a & character, no further !
!! characters are alowed to the right !
!! Function TEST_NEXT() : tests contains of next input field, return values!
!! are : E (blank); C (text); N (number); !
!! K (bracket); S (! comment); X (&END); !
!! O all other !
!! Subroutine CFIELD(string,length) : extracts the next field from the !
!! input line, at most length charcters !
!! are returned in string, if length is 0 !
!! the full field is returned, possible !
!! field seperators are: blank, comma, !
!! colon, semicolon, bracket, equal sign !
!! Subroutine P_ERROR : writes out current input line and current position !
!! within the line, does not stop code !
!! Function GET_INT() : read an integer from the input line !
!! Function GET_REAL() : read an real value from the input line !
!!<Function GET_ARRAY(I,J) : read an array field [I:J] or [I] !
!-----------------------------------------------------------------------------!
/*----------------------------------------------------------------------------|
| CP2K: A general program to perform molecular dynamics simulations |
| Copyright (C) 2024 CP2K developers group |
|-----------------------------------------------------------------------------|
|-----------------------------------------------------------------------------|
| > F R E E - F O R M A T I N P U T R O U T I N E S |
| |
| Function INSCAN(iunit,label) : looks for an input dection started by |
| label un file associated with unit iunit |
| Subroutine READ_LINE(iunit) : reads an input line, physical end of lines |
| can be masked by a & character, no further |
| characters are alowed to the right |
| Function TEST_NEXT() : tests contains of next input field, return values|
| are : E (blank); C (text); N (number); |
| K (bracket); S (! comment); X (&END); |
| O all other |
| Subroutine CFIELD(string,length) : extracts the next field from the |
| input line, at most length charcters |
| are returned in string, if length is 0 |
| the full field is returned, possible |
| field seperators are: blank, comma, |
| colon, semicolon, bracket, equal sign |
| Subroutine P_ERROR : writes out current input line and current position |
| within the line, does not stop code |
| Function GET_INT() : read an integer from the input line |
| Function GET_REAL() : read an real value from the input line |
| <Function GET_ARRAY(I,J) : read an array field [I:J] or [I] |
|----------------------------------------------------------------------------*/
MODULE parser
USE kinds, ONLY : dbl
USE global_types, ONLY : global_environment_type
USE mathconstants, ONLY : zero, one
USE stop_program, ONLY : stop_prg
USE string_utilities, ONLY : uppercase, compress
USE util, ONLY : get_unit
USE mp, ONLY : mp_bcast, mp_sync
PRIVATE
PUBLIC :: parser_init, parser_end, read_line, test_next, cfield, &
p_error, get_int, get_real, stop_parser
INTEGER, PARAMETER :: max_line = 20
CHARACTER ( LEN = 132 ) :: inp_line(max_line)
INTEGER :: num_inp_line
INTEGER :: p_line, p_char
INTEGER :: iunit
CHARACTER ( LEN = 1 ), PARAMETER :: newline = '&'
CHARACTER ( LEN = 1 ), PARAMETER :: blank = ' '
TYPE(global_environment_type) :: parenv
LOGICAL :: parallel_io, ionode
#include "global_types.h"
#include "parser.h"
#include "util.h"
CONTAINS
using namespace CP2K_NS;
!******************************************************************************
global_environment_type parenv;
SUBROUTINE parser_init ( file_name, label, ierr, globenv )
IMPLICIT NONE
INTEGER, INTENT ( OUT ) :: ierr
CHARACTER ( LEN = * ), INTENT ( IN ) :: file_name
CHARACTER ( LEN = * ), INTENT ( INOUT ) :: label
TYPE(global_environment_type), INTENT ( IN ), OPTIONAL :: globenv
CHARACTER ( LEN = 132 ) :: line
INTEGER :: ios
LOGICAL :: exists
!------------------------------------------------------------------------------
ierr = 0
! setup parameters for parallel parser
IF ( Present(globenv) ) THEN
ionode = globenv%ionode
parallel_io = .true.
parenv = globenv
ELSE
ionode = .true.
parallel_io = .false.
ENDIF
IF(ionode) THEN
!..find a not associated unit
iunit = get_unit()
line = blank
!..check the file name
INQUIRE ( FILE = file_name, EXIST = exists )
IF ( .NOT. exists ) ierr = 1
CALL uppercase ( label )
IF ( ierr == 0 ) THEN
!..look for the label
OPEN (unit=iunit,file=file_name)
REWIND (iunit)
DO
READ (iunit,iostat=ios,fmt='(a )' ) line
IF ( ios >0) THEN
CALL stop_parser ( 'parser_init','read error')
ELSE IF ( ios <0) THEN
ierr = 2
EXIT
ELSE
CALL uppercase(line)
IF (index(line,label) /= 0 ) EXIT
END IF
END DO
END IF
END IF
IF ( parallel_io ) call mp_bcast(ierr,parenv%source,parenv%group)
END SUBROUTINE parser_init
/******************************************************************************/
void parser_init(std::string file_name, std::string label, FILE ierr,
global_environment_type *globenv) {
std::string line;
FILE ios;
bool exists;
// setup parameters for parallel parser
if (globenv) {
ionode = globenv->ionode;
parallel_io = true;
parenv = *globenv;
} else {
ionode = true;
parallel_io = false;
}
if (ionode) {
//..find a not associated unit
iunit = get_unit();
line = blank;
}
}
!******************************************************************************

24
src/parser.h Normal file
View file

@ -0,0 +1,24 @@
#ifndef _PARSER_H
#define _PARSER_H
#ifdef __cplusplus
namespace CP2K_NS {
#include <string>
#endif
const int max_line = 20;
std::string inp_line[max_line];
int num_inp_line;
int p_line, p_har;
int iunit;
char newline = '&';
char blank = ' ';
bool parallel_io, ionode;
#ifdef __cplusplus
}
#endif
#endif

36
src/particle.h Normal file
View file

@ -0,0 +1,36 @@
#ifndef _PARTICLE_H
#define _PARTILCE_H
#include "linklist.h"
#ifdef __cplusplus
namespace CP2K_NS {
#endif
typedef struct particle_prop_type {
double charge;
double mass;
int ptype;
} particle_prop_type;
typedef struct particle_type {
int iatom;
double r[3];
double f[3];
double v[3];
particle_prop_type *prop;
cell_neighborlist_type *cell;
} particle_type;
typedef struct particle_list_type {
int particle_index;
double *r;
double *charge;
particle_list_type *next;
} particle_list_type;
#ifdef __cplusplus
}
#endif
#endif

0
src/particle_lists.F → src/particle_lists.c
View file

0
src/particle_types.F → src/particle_types.c
View file

0
src/periodic_table.F → src/periodic_table.c
View file

0
src/physcon.F → src/physcon.c
View file

0
src/pme.F → src/pme.c
View file

79
src/pw.h Normal file
View file

@ -0,0 +1,79 @@
#ifndef _PW_H
#define _PW_H
#ifdef __cplusplus
namespace CP2K_NS {
#else
#include <stdbool.h>
#endif
const int REALDATA1D = 301;
const int COMPLEXDATA1D = 302;
const int REALDATA3D = 303;
const int COMPLEXDATA3D = 304;
const int NOSPACE = 371;
const int REALSPACE = 372;
const int RECIPROCALSPACE = 373;
const int HALFSPACE = 211;
const int FULLSPACE = 212;
typedef struct map_pn {
int *pos;
int *neg;
} map_pn;
typedef struct pw_grid_type {
int ngpts; // # grid points
int ngpts_cut; // # grid points within cutoff
int ndimension; // # of dimensions
int bounds[3][2]; // lower and upper bounds
int npts[3]; // # point in all directions
double dr[3]; // grid spacing
double dvol, vol; // volume element, volume
struct map_pn mapl, mapm, mapn; // mapping 1D => 3D
double **g; // grid point vectors
double *gsq; // squared vector lengths
int **g_hat; // grid point indices (Miller)
// int grid_kind // type of grid: equidistance,
// logarithmic, ...
int grid_span; // type of equidistance grid:
// is f^*(g) = f(-g)
bool have_g0; // whether I have G = [0,0,0]
int first_gne0; // first g index /= 0 [1/2]
int nglengths; // # different lengths/shells
double *glsq; // squared vector lengths
int *glidx; // index to the full arrays
} pw_grid_type;
typedef struct pw_type {
double *cr;
double ***cr3d;
double *cc;
double ***cc3d;
pw_grid_type pw_grid;
int in_use; // Which data is used [r1d/c1d/r3d/c3d]
int in_space; // Real/Reciprocal space
int group_id; // In which parallel group am I
bool group_head; // Whether I am the head of this group
int group_head_id; // The group id for the head of this group
} pw_type;
typedef struct pw_grid_types {
int number;
} pw_grid_types;
void pw_zero(pw_type);
void pw_allocate(pw_type, pw_grid_type, int);
void pw_deallocate(pw_type);
#ifdef __cplusplus
}
#endif
#endif

0
src/pw_grid_types.F → src/pw_grid_types.c
View file

0
src/pw_grids.F → src/pw_grids.c
View file

175
src/pw_types.F
View file

@ -1,175 +0,0 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 CP2K developers group !
!-----------------------------------------------------------------------------!
MODULE pw_types
USE kinds, ONLY : dbl
USE pw_grid_types, ONLY : pw_grid_type
USE stop_program, ONLY : stop_prg, stop_memory
IMPLICIT NONE
PRIVATE
PUBLIC :: pw_type, pw_allocate, pw_deallocate, pw_zero
PUBLIC :: REALDATA1D, COMPLEXDATA1D, REALDATA3D, COMPLEXDATA3D
PUBLIC :: NOSPACE, REALSPACE, RECIPROCALSPACE
TYPE pw_type
REAL ( dbl ), DIMENSION ( : ), POINTER :: cr
REAL ( dbl ), DIMENSION ( :, :, : ), POINTER :: cr3d
COMPLEX ( dbl ), DIMENSION ( : ), POINTER :: cc
COMPLEX ( dbl ), DIMENSION ( :, :, : ), POINTER :: cc3d
TYPE ( pw_grid_type ), POINTER :: pw_grid
INTEGER :: in_use ! Which data is used [r1d/c1d/r3d/c3d]
INTEGER :: in_space ! Real/Reciprocal space
INTEGER :: group_id ! In which parallel group am I
LOGICAL :: group_head ! Whether I am the head of this group
INTEGER :: group_head_id ! The group id for the head of this group
END TYPE pw_type
!! Flags for the structure member 'in_use'
INTEGER, PARAMETER :: REALDATA1D = 301, COMPLEXDATA1D = 302
INTEGER, PARAMETER :: REALDATA3D = 303, COMPLEXDATA3D = 304
!! Flags for the structure member 'in_space'
INTEGER, PARAMETER :: NOSPACE = 371, REALSPACE = 372, RECIPROCALSPACE = 373
CONTAINS
!******************************************************************************
SUBROUTINE pw_allocate ( pw, pw_grid, use_data )
IMPLICIT NONE
! Arguments
INTEGER, INTENT ( IN ) :: use_data
TYPE ( pw_grid_type ), INTENT ( IN ), TARGET :: pw_grid
TYPE ( pw_type ), INTENT ( INOUT ), TARGET :: pw
! Locals
INTEGER :: allocstat
INTEGER, DIMENSION ( :, : ), POINTER :: bounds
!------------------------------------------------------------------------------
pw % in_use = use_data
pw % pw_grid => pw_grid
pw % in_space = NOSPACE
bounds => pw % pw_grid % bounds
NULLIFY ( pw % cr, pw % cc, pw % cr3d, pw % cc3d )
IF ( use_data == REALDATA1D ) THEN
ALLOCATE ( pw % cr ( pw % pw_grid % ngpts_cut ), STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "pw_allocate", "cr", pw % pw_grid % ngpts_cut )
END IF
ELSE IF ( use_data == COMPLEXDATA1D ) THEN
ALLOCATE ( pw % cc ( pw % pw_grid % ngpts_cut ), STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "pw_allocate", "cc", pw % pw_grid % ngpts_cut )
END IF
ELSE IF ( use_data == REALDATA3D ) THEN
ALLOCATE ( pw % cr3d ( &
bounds ( 1, 1 ) : bounds ( 2, 1 ), &
bounds ( 1, 2 ) : bounds ( 2, 2 ), &
bounds ( 1, 3 ) : bounds ( 2, 3 ) ), STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "pw_allocate", "cr3d", 0 )
END IF
ELSE IF ( use_data == COMPLEXDATA3D ) THEN
ALLOCATE ( pw % cc3d ( &
bounds ( 1, 1 ) : bounds ( 2, 1 ), &
bounds ( 1, 2 ) : bounds ( 2, 2 ), &
bounds ( 1, 3 ) : bounds ( 2, 3 ) ), STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "pw_allocate", "cc3d", 0 )
END IF
ELSE
CALL stop_prg ( "pw_allocate", "no possible data field" )
END IF
END SUBROUTINE pw_allocate
!******************************************************************************
SUBROUTINE pw_deallocate ( pw )
IMPLICIT NONE
! Arguments
TYPE ( pw_type ), INTENT ( INOUT ) :: pw
! Locals
INTEGER :: allocstat
!------------------------------------------------------------------------------
IF ( pw % in_use == REALDATA1D ) THEN
DEALLOCATE ( pw % cr, STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "pw_deallocate", "cr" )
END IF
ELSE IF ( pw % in_use == COMPLEXDATA1D ) THEN
DEALLOCATE ( pw % cc, STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "pw_deallocate", "cc" )
END IF
ELSE IF ( pw % in_use == REALDATA3D ) THEN
DEALLOCATE ( pw % cr3d, STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "pw_deallocate", "cr3d" )
END IF
ELSE IF ( pw % in_use == COMPLEXDATA3D ) THEN
DEALLOCATE ( pw % cc3d, STAT = allocstat )
IF ( allocstat /= 0 ) THEN
CALL stop_memory ( "pw_deallocate", "cc3d" )
END IF
ELSE
CALL stop_prg ( "pw_deallocate", "no possible data field" )
END IF
END SUBROUTINE pw_deallocate
!******************************************************************************
SUBROUTINE pw_zero ( pw )
IMPLICIT NONE
! Arguments
TYPE ( pw_type ), INTENT ( INOUT ) :: pw
!------------------------------------------------------------------------------
IF ( pw % in_use == REALDATA1D ) THEN
pw % cr = 0.0_dbl
ELSE IF ( pw % in_use == COMPLEXDATA1D ) THEN
pw % cc = 0.0_dbl
ELSE IF ( pw % in_use == REALDATA3D ) THEN
pw % cr3d = 0.0_dbl
ELSE IF ( pw % in_use == COMPLEXDATA3D ) THEN
pw % cc3d = 0.0_dbl
ELSE
CALL stop_prg ( "pw_zero", "no possible data field" )
END IF
END SUBROUTINE pw_zero
!******************************************************************************
!******************************************************************************
END MODULE pw_types

Some files were not shown because too many files have changed in this diff Show more