-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathForce_model_DoubleWell.F90
95 lines (83 loc) · 2.85 KB
/
Force_model_DoubleWell.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
! Obtaining force(fx) from position(x)
! ++++++++++++++++++++++++++++
! ++++++ Double well 1D ++++++
! ++++++++++++++++++++++++++++
subroutine Force_model_DoubleWell1D
Use Parameters, &
only: r, fr, Natom, Nbead, Eenergy, potential, &
alabel, dp_inv, address, istepsv, &
Lsave_force, &
AUtoAng, KtoAU, AngtoAU
implicit none
integer :: i, j, imode, iatom
real(8) :: f_two(3), power
integer :: Udis, Ucoor, Ufor, Uene
real(8), parameter :: xwid = 0.5 * AngtoAU
real(8), parameter :: height = 1000 * KtoAU
! +++ Calculating Forcde which atom (i) feels from atom (j) +++
fr(:,:,:) = 0.0d0
do imode = 1, Nbead
do i = 1, Natom
fr(1,i,imode) = - 4.0 * height / xwid**4 * r(1,i,imode) * (r(1,i,imode)**2 - xwid**2)
end do
end do
fr(:,:,:) = fr(:,:,:) * dp_inv
! +++ End Calculating Forcde which atom (i) feels from atom (j) +++
! +++ Calculating enetemp +++
Eenergy(:) = 0.0d0
do imode = 1, Nbead
do i = 1, Natom
Eenergy(imode) = Eenergy(imode) + height / xwid**4 * (r(1,i,imode)**2 - xwid**2)**2
end do
end do
! +++ End Calculating enetemp +++
9998 format(3E23.15)
9999 format(a2,1x,E15.9,1x,E15.9,1x,E15.9)
return
end subroutine Force_model_DoubleWell1D
! ++++++++++++++++++++++++++++
! ++++++ Double well 1D ++++++
! ++++++++++++++++++++++++++++
! ++++++++++++++++++++++++++++
! ++++++ Double well 3D ++++++
! ++++++++++++++++++++++++++++
subroutine Force_model_DoubleWell3D
Use Parameters, &
only: r, fr, Natom, Nbead, Eenergy, potential, &
alabel, dp_inv, address, istepsv, &
Lsave_force, &
AUtoAng, KtoAU, AngtoAU
implicit none
integer :: i, j, imode, iatom
real(8) :: f_two(3), power
real(8) :: rij(3), dij2, dij1
integer :: Udis, Ucoor, Ufor, Uene
real(8), parameter :: xwid = 1.0 * AngtoAU
real(8), parameter :: height = 0.003
! +++ Calculating Forcde which atom (i) feels from atom (j) +++
fr(:,:,:) = 0.0d0
Eenergy(:) = 0.0d0
do imode = 1, Nbead
do i = 1, Natom
do j = i+1, Natom
rij(:) = r(:,i,imode) - r(:,j,imode)
dij2 = dot_product(rij(:),rij(:))
dij1 = dsqrt(dij2)
! dij1 = dsqrt(dij2) - 1.0*AngtoAU ! Shift the potential
fr(:,i,imode) = - 4.0 * height / xwid**4 * (dij2 - xwid**2) * rij(:)
! f(:,i,imode) = - 4.0 * height / xwid**4 * (dij1**2 - xwid**2) * rij(:)
fr(:,j,imode) = (-1) * fr(:,i,imode)
! +++ Calculating enetemp +++
Eenergy(imode) = Eenergy(imode) + height / xwid**4 * (dij2 - xwid**2)**2
end do
end do
end do
fr(:,:,:) = fr(:,:,:) * dp_inv
! +++ End Calculating Forcde which atom (i) feels from atom (j) +++
9998 format(3E23.15)
9999 format(a2,1x,E15.9,1x,E15.9,1x,E15.9)
return
end subroutine Force_model_DoubleWell3D
! ++++++++++++++++++++++++++++
! ++++++ Double well 3D ++++++
! ++++++++++++++++++++++++++++