- #1
movezig
- 4
- 0
Hello All,
This is my first post on any forum so I will try not to get too fancy. I have been writing in fortran for about two months now and it is my first programming experience. First, I wrote a classical dynamics program for a particle moving toward a potential barrier. Then, I wrote a dynamics program to model a harmonic oscillator with a constant frequency. Then I put the two together (but uncoupled) so that there is a harmonic oscillator moving as well as a particle approaching the potential barrier in the same simulation.
Now I am trying to write a program that couples the oscillator to the motion of the particle. So basically it would be a diatomic molecule approaching some potential barrier and if it has enough energy, products will be formed (change in oscillator frequency). If there isn't enough energy in the initial conditions, the diatomic will bounce off the barrier back to reactants.
The problem: I am having trouble getting energy to be conserved in the coupled program and I think it stems in my description of the oscillator and its frequency or its derivative, but I've been through it and am stumped.
With all that said, the way I couple the oscillator to the reaction path is through a switching function, called w_q, which is coupled by transferring the oscillator motion into translational motion along the reaction path. This function is also in charge of determining the frequency of the oscillator. If reactants, the frequency is one number, and if the diatomic reacts to products, the vibrational frequency will switch to another number. All of this is dependent on the reaction path, s.
To define some terms: s = reaction path, p_s = momentum of particle,
q = oscillator displacement, p_q = oscillator momentum, H = total energy
The integration method is a 4th order Runge-Kutta and the subroutine is just the derivatives to integrate Newton's equations of motion.
My input file (cpld1_in):
1 !ntraj
100 !nprint
-6.01 !sreflect boundary condition
4.0 !sreact boundary condition
2 !nee number of energies to propagate
1.0 !alpha
1.0 !A barrier height
16.0 !m molecule mass
1.0 !mu reduced mass
0.8 !emin minimum energy attempt
1.10 !emax maximum energy attempt
-6.0 !sinit reactant position
2.0 !q_max maximum displacement
0.005 !dt
My fortran code (cpld20.f):
-------------------------------------------
program coupled20
implicit none
common / block1 / A,alpha,m,mu,w_q,q_max,phi,s_0,tutoev
real(8) :: p_q, pq_2, pq_3, pq_4
real(8) :: p_q_dot, dpqdt2, dpqdt3, dpqdt4
real(8) :: q, q_2, q_3, q_4, q_dot, dqdt2, dqdt3, dqdt4
real(8) :: alpha, A, s_0, m, mu, dt, tutoev
real(8) :: rand, sinit, q_max, sref, sreact
real(8) :: einit, emin, emax, V, K, H, w_q, phi, pi, hbar
real(8) :: sk1, sk2, sk3, sk4, psk1, psk2, psk3, psk4
real(8) :: qk1, qk2, qk3, qk4, pqk1, pqk2, pqk3, pqk4
real(8) :: s, s_2, s_3, s_4, s_dot, dsdt2, dsdt3, dsdt4
real(8) :: p_s,ps_2,ps_3,ps_4,p_s_dot,dpsdt2,dpsdt3,dpsdt4
integer :: i,ntraj,nen,nee,nstep,nprint
open(1,file='cpld1_in',status='unknown')
read(1,*)ntraj
read(1,*)nprint
read(1,*)sref
read(1,*)sreact
read(1,*)nee
read(1,*)alpha
read(1,*)A
read(1,*)m
read(1,*)mu
read(1,*)emin
read(1,*)emax
read(1,*)sinit
read(1,*)q_max
read(1,*)dt
tutoev = 1.0364314
pi = acos(-1e0)
s_0 = 0.0
open(unit=2, file='cpld20_out', action='write', status='replace')
write(2,*) 'Time(fs) s(A) Displace(A) PE(eV) KE(eV) H(eV)'
do nen=1,nee
einit=emin+(emax-emin)*float(nen-1)/float(nee-1)
call random_number(rand)
phi = 2.0*pi*rand
do i=1,ntraj
s = sinit
w_q = 0.19/(exp(alpha*(s-s_0)) + 1.0) + 0.18
p_s = sqrt(2.0*m*einit/tutoev)
q = q_max*cos(phi)
p_q = -mu*w_q*q_max*sin(phi)
nstep = 0
do while(s > sref .and. s < sreact)
if(mod(nstep,nprint) == 0) then
V = A*exp(-alpha*(s-s_0)**2) + mu*w_q*w_q*q*q/2.0
K = tutoev*p_s*p_s/(2.0*m) + p_q*p_q/(2.0*mu)
H = V + K
write(2,130)10.0*dt*nstep,s,q,V,K,H
end if
call cpld_der(s,p_s,q,p_q,s_dot,p_s_dot,q_dot,p_q_dot)
sk1 = dt*s_dot
psk1 = dt*p_s_dot
qk1 = dt*q_dot
pqk1 = dt*p_q_dot
s_2 = s + sk1/2.0
ps_2 = p_s + psk1/2.0
q_2 = q + qk1/2.0
pq_2 = p_q + pqk1/2.0
call cpld_der(s_2,ps_2,q_2,pq_2,dsdt2,dpsdt2,dqdt2,dpqdt2)
sk2 = dt*dsdt2
psk2 = dt*dpsdt2
qk2 = dt*dqdt2
pqk2 = dt*dpqdt2
s_3 = s + sk2/2.0
ps_3 = p_s + psk2/2.0
q_3 = q + qk2/2.0
pq_3 = p_q + pqk2/2.0
call cpld_der(s_3,ps_3,q_3,pq_3,dsdt3,dpsdt3,dqdt3,dpqdt3)
sk3 = dt*dsdt3
psk3 = dt*dpsdt3
qk3 = dt*dqdt3
pqk3 = dt*dpqdt3
s_4 = s + sk3
ps_4 = p_s + psk3
q_4 = q + qk3
pq_4 = p_q + pqk3
call cpld_der(s_4,ps_4,q_4,pq_4,dsdt4,dpsdt4,dqdt4,dpqdt4)
sk4 = dt*dsdt4
psk4 = dt*dpsdt4
qk4 = dt*dqdt4
pqk4 = dt*dpqdt4
s = s + (sk1 + 2.0*sk2 + 2.0*sk3 + sk4)/6.0
p_s = p_s + (psk1 + 2.0*psk2 + 2.0*psk3 + psk4)/6.0
q = q + (qk1 + 2.0*qk2 + 2.0*qk3 + qk4)/6.0
p_q = p_q + (pqk1 + 2.0*pqk2 + 2.0*pqk3 + pqk4)/6.0
w_q = 0.19/(exp(alpha*(s-s_0)) + 1.0) + 0.18
nstep = nstep + 1
end do
end do
end do
130 format(6f9.5)
close(2)
end program coupled20
subroutine cpld_der(s,ps,q,pq,sdot,psdot,qdot,pqdot)
implicit none
common / block1 / A,alpha,m,mu,w_q,q_max,phi,s_0,tutoev
real(8) :: A, alpha, m, mu, w_q, q_max, phi, s_0, tutoev
real(8) :: s, ps, q, pq, sdot, psdot, qdot, pqdot
real(8) :: dwds1,dwds2,dw2ds,psdt1,psdt2,psdt3
sdot = ps/m
dwds1 = -0.19*alpha*exp(alpha*(s-s_0))
dwds2 = dwds1/(1.0+exp(alpha*(s-s_0)))**2
dw2ds = 2.0*dwds2*w_q
psdt1 = -pq*q_max*sin(phi)*dwds2
psdt2 = -2.0*alpha*(s-s_0)*A*exp(-alpha*(s-s_0)**2)
psdt3 = 0.5*mu*q*q*dw2ds
psdot = -1.0*(psdt1+psdt2+psdt3)
qdot = pq/mu
pqdot = -mu*w_q**2*q
end subroutine cpld_der
------------------------------------------------------------------
Questions are welcome and I really appreciate the effort of anyone who has read this far.
If I should reformat my post, let me know.
Thank you!
This is my first post on any forum so I will try not to get too fancy. I have been writing in fortran for about two months now and it is my first programming experience. First, I wrote a classical dynamics program for a particle moving toward a potential barrier. Then, I wrote a dynamics program to model a harmonic oscillator with a constant frequency. Then I put the two together (but uncoupled) so that there is a harmonic oscillator moving as well as a particle approaching the potential barrier in the same simulation.
Now I am trying to write a program that couples the oscillator to the motion of the particle. So basically it would be a diatomic molecule approaching some potential barrier and if it has enough energy, products will be formed (change in oscillator frequency). If there isn't enough energy in the initial conditions, the diatomic will bounce off the barrier back to reactants.
The problem: I am having trouble getting energy to be conserved in the coupled program and I think it stems in my description of the oscillator and its frequency or its derivative, but I've been through it and am stumped.
With all that said, the way I couple the oscillator to the reaction path is through a switching function, called w_q, which is coupled by transferring the oscillator motion into translational motion along the reaction path. This function is also in charge of determining the frequency of the oscillator. If reactants, the frequency is one number, and if the diatomic reacts to products, the vibrational frequency will switch to another number. All of this is dependent on the reaction path, s.
To define some terms: s = reaction path, p_s = momentum of particle,
q = oscillator displacement, p_q = oscillator momentum, H = total energy
The integration method is a 4th order Runge-Kutta and the subroutine is just the derivatives to integrate Newton's equations of motion.
My input file (cpld1_in):
1 !ntraj
100 !nprint
-6.01 !sreflect boundary condition
4.0 !sreact boundary condition
2 !nee number of energies to propagate
1.0 !alpha
1.0 !A barrier height
16.0 !m molecule mass
1.0 !mu reduced mass
0.8 !emin minimum energy attempt
1.10 !emax maximum energy attempt
-6.0 !sinit reactant position
2.0 !q_max maximum displacement
0.005 !dt
My fortran code (cpld20.f):
-------------------------------------------
program coupled20
implicit none
common / block1 / A,alpha,m,mu,w_q,q_max,phi,s_0,tutoev
real(8) :: p_q, pq_2, pq_3, pq_4
real(8) :: p_q_dot, dpqdt2, dpqdt3, dpqdt4
real(8) :: q, q_2, q_3, q_4, q_dot, dqdt2, dqdt3, dqdt4
real(8) :: alpha, A, s_0, m, mu, dt, tutoev
real(8) :: rand, sinit, q_max, sref, sreact
real(8) :: einit, emin, emax, V, K, H, w_q, phi, pi, hbar
real(8) :: sk1, sk2, sk3, sk4, psk1, psk2, psk3, psk4
real(8) :: qk1, qk2, qk3, qk4, pqk1, pqk2, pqk3, pqk4
real(8) :: s, s_2, s_3, s_4, s_dot, dsdt2, dsdt3, dsdt4
real(8) :: p_s,ps_2,ps_3,ps_4,p_s_dot,dpsdt2,dpsdt3,dpsdt4
integer :: i,ntraj,nen,nee,nstep,nprint
open(1,file='cpld1_in',status='unknown')
read(1,*)ntraj
read(1,*)nprint
read(1,*)sref
read(1,*)sreact
read(1,*)nee
read(1,*)alpha
read(1,*)A
read(1,*)m
read(1,*)mu
read(1,*)emin
read(1,*)emax
read(1,*)sinit
read(1,*)q_max
read(1,*)dt
tutoev = 1.0364314
pi = acos(-1e0)
s_0 = 0.0
open(unit=2, file='cpld20_out', action='write', status='replace')
write(2,*) 'Time(fs) s(A) Displace(A) PE(eV) KE(eV) H(eV)'
do nen=1,nee
einit=emin+(emax-emin)*float(nen-1)/float(nee-1)
call random_number(rand)
phi = 2.0*pi*rand
do i=1,ntraj
s = sinit
w_q = 0.19/(exp(alpha*(s-s_0)) + 1.0) + 0.18
p_s = sqrt(2.0*m*einit/tutoev)
q = q_max*cos(phi)
p_q = -mu*w_q*q_max*sin(phi)
nstep = 0
do while(s > sref .and. s < sreact)
if(mod(nstep,nprint) == 0) then
V = A*exp(-alpha*(s-s_0)**2) + mu*w_q*w_q*q*q/2.0
K = tutoev*p_s*p_s/(2.0*m) + p_q*p_q/(2.0*mu)
H = V + K
write(2,130)10.0*dt*nstep,s,q,V,K,H
end if
call cpld_der(s,p_s,q,p_q,s_dot,p_s_dot,q_dot,p_q_dot)
sk1 = dt*s_dot
psk1 = dt*p_s_dot
qk1 = dt*q_dot
pqk1 = dt*p_q_dot
s_2 = s + sk1/2.0
ps_2 = p_s + psk1/2.0
q_2 = q + qk1/2.0
pq_2 = p_q + pqk1/2.0
call cpld_der(s_2,ps_2,q_2,pq_2,dsdt2,dpsdt2,dqdt2,dpqdt2)
sk2 = dt*dsdt2
psk2 = dt*dpsdt2
qk2 = dt*dqdt2
pqk2 = dt*dpqdt2
s_3 = s + sk2/2.0
ps_3 = p_s + psk2/2.0
q_3 = q + qk2/2.0
pq_3 = p_q + pqk2/2.0
call cpld_der(s_3,ps_3,q_3,pq_3,dsdt3,dpsdt3,dqdt3,dpqdt3)
sk3 = dt*dsdt3
psk3 = dt*dpsdt3
qk3 = dt*dqdt3
pqk3 = dt*dpqdt3
s_4 = s + sk3
ps_4 = p_s + psk3
q_4 = q + qk3
pq_4 = p_q + pqk3
call cpld_der(s_4,ps_4,q_4,pq_4,dsdt4,dpsdt4,dqdt4,dpqdt4)
sk4 = dt*dsdt4
psk4 = dt*dpsdt4
qk4 = dt*dqdt4
pqk4 = dt*dpqdt4
s = s + (sk1 + 2.0*sk2 + 2.0*sk3 + sk4)/6.0
p_s = p_s + (psk1 + 2.0*psk2 + 2.0*psk3 + psk4)/6.0
q = q + (qk1 + 2.0*qk2 + 2.0*qk3 + qk4)/6.0
p_q = p_q + (pqk1 + 2.0*pqk2 + 2.0*pqk3 + pqk4)/6.0
w_q = 0.19/(exp(alpha*(s-s_0)) + 1.0) + 0.18
nstep = nstep + 1
end do
end do
end do
130 format(6f9.5)
close(2)
end program coupled20
subroutine cpld_der(s,ps,q,pq,sdot,psdot,qdot,pqdot)
implicit none
common / block1 / A,alpha,m,mu,w_q,q_max,phi,s_0,tutoev
real(8) :: A, alpha, m, mu, w_q, q_max, phi, s_0, tutoev
real(8) :: s, ps, q, pq, sdot, psdot, qdot, pqdot
real(8) :: dwds1,dwds2,dw2ds,psdt1,psdt2,psdt3
sdot = ps/m
dwds1 = -0.19*alpha*exp(alpha*(s-s_0))
dwds2 = dwds1/(1.0+exp(alpha*(s-s_0)))**2
dw2ds = 2.0*dwds2*w_q
psdt1 = -pq*q_max*sin(phi)*dwds2
psdt2 = -2.0*alpha*(s-s_0)*A*exp(-alpha*(s-s_0)**2)
psdt3 = 0.5*mu*q*q*dw2ds
psdot = -1.0*(psdt1+psdt2+psdt3)
qdot = pq/mu
pqdot = -mu*w_q**2*q
end subroutine cpld_der
------------------------------------------------------------------
Questions are welcome and I really appreciate the effort of anyone who has read this far.
If I should reformat my post, let me know.
Thank you!