Skip to content

Commit

Permalink
Remove obsolete auxiliary math routines
Browse files Browse the repository at this point in the history
  • Loading branch information
@mkrack
mkrack committed Jul 1, 2022
1 parent 8cd0e9d commit cc2b7bf
Show file tree
Hide file tree
Showing 23 changed files with 288 additions and 416 deletions.
78 changes: 2 additions & 76 deletions src/common/mathlib.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,15 +55,12 @@ MODULE mathlib
inv_3x3, &
lcm, &
vector_product, &
matmul_3x3, &
matvec_3x3, &
pswitch, &
rotate_vector, &
reflect_vector, &
dotprod_3d, &
transpose_3d, &
expint, abnormal_value, &
get_diag, set_diag
get_diag, &
set_diag

INTERFACE det_3x3
MODULE PROCEDURE det_3x3_1, det_3x3_2
Expand Down Expand Up @@ -1400,77 +1397,6 @@ FUNCTION ei(x)

END FUNCTION ei

! **************************************************************************************************
!> \brief ...
!> \param mat1 ...
!> \param mat2 ...
!> \return ...
! **************************************************************************************************
PURE FUNCTION matmul_3x3(mat1, mat2)
REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: mat1, mat2
REAL(KIND=dp), DIMENSION(3, 3) :: matmul_3x3

matmul_3x3(1, 1) = mat1(1, 1)*mat2(1, 1) + mat1(1, 2)*mat2(2, 1) + mat1(1, 3)*mat2(3, 1)
matmul_3x3(1, 2) = mat1(1, 1)*mat2(1, 2) + mat1(1, 2)*mat2(2, 2) + mat1(1, 3)*mat2(3, 2)
matmul_3x3(1, 3) = mat1(1, 1)*mat2(1, 3) + mat1(1, 2)*mat2(2, 3) + mat1(1, 3)*mat2(3, 3)
matmul_3x3(2, 1) = mat1(2, 1)*mat2(1, 1) + mat1(2, 2)*mat2(2, 1) + mat1(2, 3)*mat2(3, 1)
matmul_3x3(2, 2) = mat1(2, 1)*mat2(1, 2) + mat1(2, 2)*mat2(2, 2) + mat1(2, 3)*mat2(3, 2)
matmul_3x3(2, 3) = mat1(2, 1)*mat2(1, 3) + mat1(2, 2)*mat2(2, 3) + mat1(2, 3)*mat2(3, 3)
matmul_3x3(3, 1) = mat1(3, 1)*mat2(1, 1) + mat1(3, 2)*mat2(2, 1) + mat1(3, 3)*mat2(3, 1)
matmul_3x3(3, 2) = mat1(3, 1)*mat2(1, 2) + mat1(3, 2)*mat2(2, 2) + mat1(3, 3)*mat2(3, 2)
matmul_3x3(3, 3) = mat1(3, 1)*mat2(1, 3) + mat1(3, 2)*mat2(2, 3) + mat1(3, 3)*mat2(3, 3)
END FUNCTION matmul_3x3

! **************************************************************************************************
!> \brief ...
!> \param res ...
!> \param mat ...
!> \param vec ...
! **************************************************************************************************
PURE SUBROUTINE matvec_3x3(res, mat, vec)
REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: res
REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: mat
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vec

res(1) = mat(1, 1)*vec(1) + mat(1, 2)*vec(2) + mat(1, 3)*vec(3)
res(2) = mat(2, 1)*vec(1) + mat(2, 2)*vec(2) + mat(2, 3)*vec(3)
res(3) = mat(3, 1)*vec(1) + mat(3, 2)*vec(2) + mat(3, 3)*vec(3)
END SUBROUTINE matvec_3x3

! **************************************************************************************************
!> \brief ...
!> \param vec1 ...
!> \param vec2 ...
!> \return ...
! **************************************************************************************************
PURE FUNCTION dotprod_3d(vec1, vec2)
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: vec1, vec2
REAL(KIND=dp) :: dotprod_3d

dotprod_3d = &
vec1(1)*vec2(1) &
+ vec1(2)*vec2(2) &
+ vec1(3)*vec2(3)
END FUNCTION dotprod_3d

! **************************************************************************************************
!> \brief ...
!> \param mat ...
!> \return ...
! **************************************************************************************************
PURE FUNCTION transpose_3d(mat)
REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: mat
REAL(KIND=dp), DIMENSION(3, 3) :: transpose_3d

INTEGER :: i, j

DO i = 1, 3
DO j = 1, 3
transpose_3d(j, i) = mat(i, j)
END DO
END DO
END FUNCTION transpose_3d

! **************************************************************************************************
!> \brief computes the exponential integral
!> En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,1,..
Expand Down
133 changes: 61 additions & 72 deletions src/constraint_3x3.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,6 @@ MODULE constraint_3x3
USE constraint_fxd, ONLY: check_fixed_atom_cns_g3x3
USE kinds, ONLY: dp
USE linear_systems, ONLY: solve_system
USE mathlib, ONLY: dotprod_3d,&
matvec_3x3
USE molecule_kind_types, ONLY: fixd_constraint_type,&
g3x3_constraint_type,&
get_molecule_kind,&
Expand Down Expand Up @@ -407,15 +405,15 @@ SUBROUTINE shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
index_a, index_b, index_c, fixd_list, lg3x3(iconst))
! construct matrix
amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r0_12, lg3x3(iconst)%fa)
amat(1, 2) = imass1*DOTPROD_3D(r0_12, lg3x3(iconst)%fb)
amat(1, 3) = -imass2*DOTPROD_3D(r0_12, lg3x3(iconst)%fc)
amat(2, 1) = imass1*DOTPROD_3D(r0_13, lg3x3(iconst)%fa)
amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r0_13, lg3x3(iconst)%fb)
amat(2, 3) = imass3*DOTPROD_3D(r0_13, lg3x3(iconst)%fc)
amat(3, 1) = -imass2*DOTPROD_3D(r0_23, lg3x3(iconst)%fa)
amat(3, 2) = imass3*DOTPROD_3D(r0_23, lg3x3(iconst)%fb)
amat(3, 3) = (imass3 + imass2)*DOTPROD_3D(r0_23, lg3x3(iconst)%fc)
amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r0_12, lg3x3(iconst)%fa)
amat(1, 2) = imass1*DOT_PRODUCT(r0_12, lg3x3(iconst)%fb)
amat(1, 3) = -imass2*DOT_PRODUCT(r0_12, lg3x3(iconst)%fc)
amat(2, 1) = imass1*DOT_PRODUCT(r0_13, lg3x3(iconst)%fa)
amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r0_13, lg3x3(iconst)%fb)
amat(2, 3) = imass3*DOT_PRODUCT(r0_13, lg3x3(iconst)%fc)
amat(3, 1) = -imass2*DOT_PRODUCT(r0_23, lg3x3(iconst)%fa)
amat(3, 2) = imass3*DOT_PRODUCT(r0_23, lg3x3(iconst)%fb)
amat(3, 3) = (imass3 + imass2)*DOT_PRODUCT(r0_23, lg3x3(iconst)%fc)
! Store values
lg3x3(iconst)%r0_12 = r0_12
lg3x3(iconst)%r0_13 = r0_13
Expand All @@ -438,17 +436,17 @@ SUBROUTINE shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
lg3x3(iconst)%lambda(2)*imass1*lg3x3(iconst)%fb - &
lg3x3(iconst)%lambda(3)*imass2*lg3x3(iconst)%fc
bvec(1, 1) = g3x3_list(iconst)%dab*g3x3_list(iconst)%dab &
- dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_12, r0_12)
- dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_12, r0_12)
vec = lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa*imass1 + &
lg3x3(iconst)%lambda(2)*(imass1 + imass3)*lg3x3(iconst)%fb + &
lg3x3(iconst)%lambda(3)*imass3*lg3x3(iconst)%fc
bvec(2, 1) = g3x3_list(iconst)%dac*g3x3_list(iconst)%dac &
- dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_13, r0_13)
- dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_13, r0_13)
vec = -lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa*imass2 + &
lg3x3(iconst)%lambda(2)*imass3*lg3x3(iconst)%fb + &
lg3x3(iconst)%lambda(3)*(imass2 + imass3)*lg3x3(iconst)%fc
bvec(3, 1) = g3x3_list(iconst)%dbc*g3x3_list(iconst)%dbc &
- dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_23, r0_23)
- dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_23, r0_23)
bvec = bvec*idtsq
! get lambda
atemp = amat
Expand Down Expand Up @@ -566,19 +564,19 @@ SUBROUTINE shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
index_a, index_b, index_c, fixd_list, lg3x3(iconst))
! rotate fconst:
CALL matvec_3x3(f_roll1, r_shake, lg3x3(iconst)%fa)
CALL matvec_3x3(f_roll2, r_shake, lg3x3(iconst)%fb)
CALL matvec_3x3(f_roll3, r_shake, lg3x3(iconst)%fc)
f_roll1 = MATMUL(r_shake, lg3x3(iconst)%fa)
f_roll2 = MATMUL(r_shake, lg3x3(iconst)%fb)
f_roll3 = MATMUL(r_shake, lg3x3(iconst)%fc)
! construct matrix
amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r0_12, f_roll1)
amat(1, 2) = imass1*DOTPROD_3D(r0_12, f_roll2)
amat(1, 3) = -imass2*DOTPROD_3D(r0_12, f_roll3)
amat(2, 1) = imass1*DOTPROD_3D(r0_13, f_roll1)
amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r0_13, f_roll2)
amat(2, 3) = imass3*DOTPROD_3D(r0_13, f_roll3)
amat(3, 1) = -imass2*DOTPROD_3D(r0_23, f_roll1)
amat(3, 2) = imass3*DOTPROD_3D(r0_23, f_roll2)
amat(3, 3) = (imass3 + imass2)*DOTPROD_3D(r0_23, f_roll3)
amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r0_12, f_roll1)
amat(1, 2) = imass1*DOT_PRODUCT(r0_12, f_roll2)
amat(1, 3) = -imass2*DOT_PRODUCT(r0_12, f_roll3)
amat(2, 1) = imass1*DOT_PRODUCT(r0_13, f_roll1)
amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r0_13, f_roll2)
amat(2, 3) = imass3*DOT_PRODUCT(r0_13, f_roll3)
amat(3, 1) = -imass2*DOT_PRODUCT(r0_23, f_roll1)
amat(3, 2) = imass3*DOT_PRODUCT(r0_23, f_roll2)
amat(3, 3) = (imass3 + imass2)*DOT_PRODUCT(r0_23, f_roll3)
! Store values
lg3x3(iconst)%r0_12 = r0_12
lg3x3(iconst)%r0_13 = r0_13
Expand Down Expand Up @@ -607,17 +605,17 @@ SUBROUTINE shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
lg3x3(iconst)%lambda(2)*imass1*f_roll2 - &
lg3x3(iconst)%lambda(3)*imass2*f_roll3
bvec(1, 1) = g3x3_list(iconst)%dab*g3x3_list(iconst)%dab &
- dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_12, r0_12)
- dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_12, r0_12)
vec = lg3x3(iconst)%lambda(1)*f_roll1*imass1 + &
lg3x3(iconst)%lambda(2)*(imass1 + imass3)*f_roll2 + &
lg3x3(iconst)%lambda(3)*imass3*f_roll3
bvec(2, 1) = g3x3_list(iconst)%dac*g3x3_list(iconst)%dac &
- dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_13, r0_13)
- dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_13, r0_13)
vec = -lg3x3(iconst)%lambda(1)*f_roll1*imass2 + &
lg3x3(iconst)%lambda(2)*imass3*f_roll2 + &
lg3x3(iconst)%lambda(3)*(imass2 + imass3)*f_roll3
bvec(3, 1) = g3x3_list(iconst)%dbc*g3x3_list(iconst)%dbc &
- dtsqby2*dtsqby2*DOTPROD_3D(vec, vec) - DOTPROD_3D(r0_23, r0_23)
- dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_23, r0_23)
bvec = bvec*idtsq

! get lambda
Expand All @@ -634,18 +632,12 @@ SUBROUTINE shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
fc3 = -lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb - &
lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
CALL MATVEC_3x3(vec, r_shake, fc1)
r1(:) = pos(:, index_a) + imass1*dtsqby2*vec
CALL MATVEC_3x3(vec, r_shake, fc2)
r2(:) = pos(:, index_b) + imass2*dtsqby2*vec
CALL MATVEC_3x3(vec, r_shake, fc3)
r3(:) = pos(:, index_c) + imass3*dtsqby2*vec
CALL MATVEC_3x3(vec, v_shake, fc1)
v1(:) = vel(:, index_a) + imass1*dtby2*vec
CALL MATVEC_3x3(vec, v_shake, fc2)
v2(:) = vel(:, index_b) + imass2*dtby2*vec
CALL MATVEC_3x3(vec, v_shake, fc3)
v3(:) = vel(:, index_c) + imass3*dtby2*vec
r1(:) = pos(:, index_a) + imass1*dtsqby2*MATMUL(r_shake, fc1)
r2(:) = pos(:, index_b) + imass2*dtsqby2*MATMUL(r_shake, fc2)
r3(:) = pos(:, index_c) + imass3*dtsqby2*MATMUL(r_shake, fc3)
v1(:) = vel(:, index_a) + imass1*dtby2*MATMUL(v_shake, fc1)
v2(:) = vel(:, index_b) + imass2*dtby2*MATMUL(v_shake, fc2)
v3(:) = vel(:, index_c) + imass3*dtby2*MATMUL(v_shake, fc3)
r12 = r1 - r2
r13 = r1 - r3
r23 = r2 - r3
Expand Down Expand Up @@ -703,7 +695,7 @@ SUBROUTINE rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
REAL(KIND=dp) :: dtby2, idt, imass1, imass2, imass3, mass
REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, fc1, fc2, &
fc3, lambda, r12, r13, r23, v12, v13, &
v23, vec
v23
REAL(KIND=dp), DIMENSION(3, 1) :: bvec
REAL(KIND=dp), DIMENSION(3, 3) :: amat
TYPE(atomic_kind_type), POINTER :: atomic_kind
Expand Down Expand Up @@ -737,28 +729,25 @@ SUBROUTINE rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
index_a, index_b, index_c, fixd_list, lg3x3(iconst))
! roll the fc
CALL MATVEC_3x3(f_roll1, r_rattle, lg3x3(iconst)%fa)
CALL MATVEC_3x3(f_roll2, r_rattle, lg3x3(iconst)%fb)
CALL MATVEC_3x3(f_roll3, r_rattle, lg3x3(iconst)%fc)
f_roll1 = MATMUL(r_rattle, lg3x3(iconst)%fa)
f_roll2 = MATMUL(r_rattle, lg3x3(iconst)%fb)
f_roll3 = MATMUL(r_rattle, lg3x3(iconst)%fc)

! construct matrix
amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r12, f_roll1)
amat(1, 2) = imass1*DOTPROD_3D(r12, f_roll2)
amat(1, 3) = -imass2*DOTPROD_3D(r12, f_roll3)
amat(2, 1) = imass1*DOTPROD_3D(r13, f_roll1)
amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r13, f_roll2)
amat(2, 3) = imass3*DOTPROD_3D(r13, f_roll3)
amat(3, 1) = -imass2*DOTPROD_3D(r23, f_roll1)
amat(3, 2) = imass3*DOTPROD_3D(r23, f_roll2)
amat(3, 3) = (imass2 + imass3)*DOTPROD_3D(r23, f_roll3)
amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r12, f_roll1)
amat(1, 2) = imass1*DOT_PRODUCT(r12, f_roll2)
amat(1, 3) = -imass2*DOT_PRODUCT(r12, f_roll3)
amat(2, 1) = imass1*DOT_PRODUCT(r13, f_roll1)
amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r13, f_roll2)
amat(2, 3) = imass3*DOT_PRODUCT(r13, f_roll3)
amat(3, 1) = -imass2*DOT_PRODUCT(r23, f_roll1)
amat(3, 2) = imass3*DOT_PRODUCT(r23, f_roll2)
amat(3, 3) = (imass2 + imass3)*DOT_PRODUCT(r23, f_roll3)

! construct solution vector
CALL matvec_3x3(vec, veps, r12)
bvec(1, 1) = DOTPROD_3D(r12, v12 + vec)
CALL matvec_3x3(vec, veps, r13)
bvec(2, 1) = DOTPROD_3D(r13, v13 + vec)
CALL matvec_3x3(vec, veps, r23)
bvec(3, 1) = DOTPROD_3D(r23, v23 + vec)
bvec(1, 1) = DOT_PRODUCT(r12, v12 + MATMUL(veps, r12))
bvec(2, 1) = DOT_PRODUCT(r13, v13 + MATMUL(veps, r13))
bvec(3, 1) = DOT_PRODUCT(r23, v23 + MATMUL(veps, r23))
bvec = -bvec*2.0_dp*idt

! get lambda
Expand Down Expand Up @@ -838,20 +827,20 @@ SUBROUTINE rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
index_a, index_b, index_c, fixd_list, lg3x3(iconst))
! construct matrix
amat(1, 1) = (imass1 + imass2)*DOTPROD_3D(r12, lg3x3(iconst)%fa)
amat(1, 2) = imass1*DOTPROD_3D(r12, lg3x3(iconst)%fb)
amat(1, 3) = -imass2*DOTPROD_3D(r12, lg3x3(iconst)%fc)
amat(2, 1) = imass1*DOTPROD_3D(r13, lg3x3(iconst)%fa)
amat(2, 2) = (imass1 + imass3)*DOTPROD_3D(r13, lg3x3(iconst)%fb)
amat(2, 3) = imass3*DOTPROD_3D(r13, lg3x3(iconst)%fc)
amat(3, 1) = -imass2*DOTPROD_3D(r23, lg3x3(iconst)%fa)
amat(3, 2) = imass3*DOTPROD_3D(r23, lg3x3(iconst)%fb)
amat(3, 3) = (imass2 + imass3)*DOTPROD_3D(r23, lg3x3(iconst)%fc)
amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r12, lg3x3(iconst)%fa)
amat(1, 2) = imass1*DOT_PRODUCT(r12, lg3x3(iconst)%fb)
amat(1, 3) = -imass2*DOT_PRODUCT(r12, lg3x3(iconst)%fc)
amat(2, 1) = imass1*DOT_PRODUCT(r13, lg3x3(iconst)%fa)
amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r13, lg3x3(iconst)%fb)
amat(2, 3) = imass3*DOT_PRODUCT(r13, lg3x3(iconst)%fc)
amat(3, 1) = -imass2*DOT_PRODUCT(r23, lg3x3(iconst)%fa)
amat(3, 2) = imass3*DOT_PRODUCT(r23, lg3x3(iconst)%fb)
amat(3, 3) = (imass2 + imass3)*DOT_PRODUCT(r23, lg3x3(iconst)%fc)

! construct solution vector
bvec(1, 1) = DOTPROD_3D(r12, v12)
bvec(2, 1) = DOTPROD_3D(r13, v13)
bvec(3, 1) = DOTPROD_3D(r23, v23)
bvec(1, 1) = DOT_PRODUCT(r12, v12)
bvec(2, 1) = DOT_PRODUCT(r13, v13)
bvec(3, 1) = DOT_PRODUCT(r23, v23)
bvec = -bvec*2.0_dp*idt

! get lambda
Expand Down

0 comments on commit cc2b7bf

Please sign in to comment.