Module that contains some useful utilities. More...

Functions/Subroutines
real(8) function, dimension(3)	pbc (a, box)
	Corrects for the periodic boundary condition. More...

real(8) function, dimension(3)	cross (a, b)
	Performs cross product between two vectors. More...

real(8) function	distance2 (a, b, box)
	Calculates the distance squared between two points. More...

real(8) function	distance (a, b, box)
	Calculates the distance between two points. More...

real(8) function, dimension(3)	bond_vector (a, b, box)
	Calculates the bond vector between two points. More...

real(8) function	magnitude (a)
	Calculates the magnitude of a vector. More...

real(8) function	bond_angle (a, b, c, box)
	Calculates the bond angle between two vectors. More...

real(8) function	dihedral_angle (i, j, k, l, box)
	Calculates the dihedral angle between two planes formed by four atoms. More...

Detailed Description

Module that contains some useful utilities.

Function/Subroutine Documentation

◆ bond_angle()

real(8) function dcdfort_utils::bond_angle	(	real(8), dimension(3), intent(in)	a,
		real(8), dimension(3), intent(in)	b,
		real(8), dimension(3), intent(in)	c,
		real(8), dimension(6), intent(in)	box
	)

Calculates the bond angle between two vectors.

Calculates the angle between the vector formed by a-b and b-c.

Parameters

[in]	a	first point
[in]	b	middle point
[in]	c	third point
[in]	box	box, if pbc to be accounted for

Returns: bond angle between a-b-c

 
         implicit none
         real(8) :: bond_angle
         real(8), intent(in), dimension(3) :: a, b, c
         real(8), intent(in) :: box(6)
         real(8), dimension(3) :: bond1, bond2
 
         bond1 = bond_vector(b, a, box)
         bond2 = bond_vector(b, c, box)
 
         bond_angle = dacos(dot_product(bond1, bond2)/(magnitude(bond1)*magnitude(bond2)))
 

◆ bond_vector()

real(8) function, dimension(3) dcdfort_utils::bond_vector	(	real(8), dimension(3), intent(in)	a,
		real(8), dimension(3), intent(in)	b,
		real(8), dimension(6), intent(in)	box
	)

Calculates the bond vector between two points.

Parameters

[in]	a	first point
[in]	b	second point
[in]	box	box, if pbc to be accounted for

Returns: bond vector pointing from a to b

 
         implicit none
         real(8) :: bond_vector(3)
         real(8), intent(in), dimension(3) :: a, b
         real(8), intent(in) :: box(6)
 
         bond_vector = pbc(a-b, box)
 

◆ cross()

real(8) function, dimension(3) dcdfort_utils::cross	(	real(8), dimension(3), intent(in)	a,
		real(8), dimension(3), intent(in)	b
	)

Performs cross product between two vectors.

Parameters

[in]	a	first vector
[in]	b	second vector

Returns: resulting cross product

 
         implicit none
         real(8) :: cross(3)
         real(8), intent(in), dimension(3) :: a, b
 
         cross(1) = a(2) * b(3) - a(3) * b(2)
         cross(2) = a(3) * b(1) - a(1) * b(3)
         cross(3) = a(1) * b(2) - a(2) * b(1)
 

◆ dihedral_angle()

real(8) function dcdfort_utils::dihedral_angle	(	real(8), dimension(3), intent(in)	i,
		real(8), dimension(3), intent(in)	j,
		real(8), dimension(3), intent(in)	k,
		real(8), dimension(3), intent(in)	l,
		real(8), dimension(6), intent(in)	box
	)

Calculates the dihedral angle between two planes formed by four atoms.

Calculates the dihedral angle between the vectors formed by i-j, j-k, k-l

Parameters

[in]	i	first point
[in]	j	middle point
[in]	k	third point
[in]	l	fourth point
[in]	box	box, if pbc to be accounted for

Returns: dihedral angle forms by i-j-k-l

 
         implicit none
         real(8) :: dihedral_angle
         real(8), intent(in), dimension(3) :: i, j, k, l
         real(8), intent(in), dimension(6) :: box
         real(8) :: a_mag, b_mag, g_mag
         real(8), dimension(3) :: h, g, f, a, b, cross_ba
 
         h = bond_vector(k, l, box)
         g = bond_vector(k, j, box)
         f = bond_vector(j, i, box)
         a = cross(f, g)
         b = cross(h, g)
         cross_ba = cross(b, a)
         a_mag = magnitude(a)
         b_mag = magnitude(b)
         g_mag = magnitude(g)
 
         dihedral_angle = atan2(dot_product(cross_ba, g) / (a_mag*b_mag*g_mag), dot_product(a, b) / (a_mag*b_mag))
 

◆ distance()

real(8) function dcdfort_utils::distance	(	real(8), dimension(3), intent(in)	a,
		real(8), dimension(3), intent(in)	b,
		real(8), dimension(6), intent(in), optional	box
	)

Calculates the distance between two points.

Parameters

[in]	a	first point
[in]	b	second point
[in]	box	box, if pbc to be accounted for

Returns: distance

 
         implicit none
         real(8) :: distance
         real(8), intent(in), dimension(3) :: a, b
         real(8), intent(in), optional :: box(6)
 
         if (present(box)) then
             distance = dsqrt(distance2(a, b, box))
         else
             distance = dsqrt(distance2(a, b))
         end if
 

◆ distance2()

real(8) function dcdfort_utils::distance2	(	real(8), dimension(3), intent(in)	a,
		real(8), dimension(3), intent(in)	b,
		real(8), dimension(6), intent(in), optional	box
	)

Calculates the distance squared between two points.

Parameters

[in]	a	first point
[in]	b	second point
[in]	box	box, if pbc to be accounted for

Returns: distance squared

 
         implicit none
         real(8) :: distance2
         real(8), intent(in), dimension(3) :: a, b
         real(8) :: c(3)
         real(8), intent(in), optional :: box(6)
 
         if (present(box)) then
             c = pbc(a - b, box)
         else
             c = a - b
         end if
         distance2 = dot_product(c, c)
 

◆ magnitude()

real(8) function dcdfort_utils::magnitude ( real(8), dimension(3), intent(in) a )

Calculates the magnitude of a vector.

Parameters

[in] a vector

Returns: magnitude of a

 
         implicit none
         real(8) :: magnitude
         real(8), intent(in) :: a(3)
 
         magnitude = dsqrt(dot_product(a, a))

◆ pbc()

real(8) function, dimension(3) dcdfort_utils::pbc	(	real(8), dimension(3), intent(in)	a,
		real(8), dimension(6), intent(in)	box
	)

Corrects for the periodic boundary condition.

Moves particle (or vector) the distance of half the box if it is more than half the distance of the box

Parameters

[in]	a	original coordinates
[in]	box	simulation box

Returns: the shifted coordinates

 
         implicit none
         real(8), intent(in) :: a(3), box(6)
         real(8) :: pbc(3), tbox(3,3) = 0.0d0
         integer :: i, shift
 
         ! A = box(1), length of unit cell vector along x-axis;
         ! B = box(2), length of unit cell vector in xy-plane;
         ! C = box(3), length of unit cell vector in yz-plane;
         ! alpha = box(4), cosine of angle between B and C
         ! beta  = box(5), cosine of angle between A and C
         ! gamma = box(6), cosine of angle between A and B
 
         ! New matrix made up of box vectors:
 
         !  ax  bx  cx
         !   0  by  cy
         !   0   0  cz
 
         ! convert angles to box vectors
         
         ! A**2 = ax**2 + ay**2 + az**2
         ! A**2 = ax**2 + (0)**2 + (0)**2
         ! ax = A
         tbox(1,1) = box(1)
 
         ! dot(A_vec,B_vec) = ax*bx + ay*by + az*bz
         ! dot(A_vec,B_vec) = ax*bx + (0)*by + (0)*bz
         ! A*B*gamma = ax*bx
         ! bx = B*gamma
         tbox(1,2) = box(2)*box(6)
 
         ! dot(A_vec,C_vec) = ax*cx + ay*cy + az*cz
         ! dot(A_vec,C_vec) = ax*cx + (0)*cy + (0)*cz
         ! A*C*beta = ax*cx
         ! cx = C*beta
         tbox(1,3) = box(3)*box(5)
 
         ! B**2 = bx**2 + by**2 + bz**2
         ! B**2 = bx**2 + by**2 + (0)**2
         ! by = sqrt(B**2 - bx**2)
         tbox(2,2) = dsqrt(box(2)**2-tbox(1,2)**2)
 
         ! dot(B_vec,C_vec) = bx*cx + by*cy + (0)*cz
         ! B*C*alpha = bx*cx + by*cy
         ! cy = (B*C*alpha - bx*cx) / by
         tbox(2,3) = (box(2)*box(3)*box(4) - tbox(1,2)*tbox(1,3))/tbox(2,2)
 
         ! C = cx**2 + cy**2 + cz**2
         ! cz = sqrt(C**2 - cx**2 - cy**2)
         tbox(3,3) = dsqrt(box(3)**2-tbox(1,3)**2-tbox(2,3)**2)
 
         pbc = a
 
         ! Now perform PBC calculation
         shift = nint(pbc(3) / tbox(3,3))
         if (shift .ne. 0) then
             pbc(3) = pbc(3) - tbox(3,3) * shift
             pbc(2) = pbc(2) - tbox(2,3) * shift
             pbc(1) = pbc(1) - tbox(1,3) * shift
         end if
 
         shift = nint(pbc(2) / tbox(2,2))
         if (shift .ne. 0) then
             pbc(2) = pbc(2) - tbox(2,2) * shift
             pbc(1) = pbc(1) - tbox(1,2) * shift
         end if
 
         shift = nint(pbc(1) / tbox(1,1))
         if (shift .ne. 0) then
             pbc(1) = pbc(1) - tbox(1,1) * shift
         end if
 

Functions/Subroutines

Detailed Description

Function/Subroutine Documentation

◆ bond_angle()

◆ bond_vector()

◆ cross()

◆ dihedral_angle()

◆ distance()

◆ distance2()

◆ magnitude()

◆ pbc()