Hello All,(adsbygoogle = window.adsbygoogle || []).push({});

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!

**Physics Forums | Science Articles, Homework Help, Discussion**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Fortran quassiclassical trajectory help

**Physics Forums | Science Articles, Homework Help, Discussion**