time = 0.0d0
do
if (time .ge. time_end) exit
if (time <= time_center) then
e_rabi_0 = e_rabi_i*exp(-(time-time_center)**2/time_width**2) ! = exp(-4*(time-time_center)**2/time_center**2)
else
e_rabi_0 = e_rabi_i
end if
e_rabi_1 = e_rabi_0
call exci_rates(de,sq_fcf,gam,kop,e_rabi_0,e_rabi_1)
do mc = 0, 1
do me = 0, 1
do mv = 0, mv_max
do nc = 0, 1
do ne = 0, 1
do nv = 0, mv_max
k_ref(mc,me,mv,nc,ne,nv) = kivr(mc,me,mv,nc,ne,nv) + kop(mc,me,mv,nc,ne,nv) + klu(mc,me,mv,nc,ne,nv)
end do
end do
end do
end do
end do
end do
voltage = voltage_give
call mol_lead_rates(voltage,de,v_l,v_r,sq_fcf,klw,klw_l)
do mc = 0, 1
do me = 0, 1
do mv = 0, mv_max
do nc = 0, 1
do ne = 0, 1
do nv = 0, mv_max
k_tot(mc,me,mv,nc,ne,nv) = k_ref(mc,me,mv,nc,ne,nv) + klw(mc,me,mv,nc,ne,nv)
end do
end do
end do
end do
end do
end do
if (time == 0.0d0) then
do mc = 0, 1
do me = 0, 1
do mv = 0, mv_max
pop(mc,me,mv) = 0.d0
end do
end do
end do
pop(0,0,0) = 1.d0
stop_time = time_final
100 do mc = 0, 1
do me = 0, 1
do mv = 0, mv_max
pop_old(mc,me,mv) = pop(mc,me,mv)
end do
end do
end do
call multi_to_single(pop_old,fin)
call rk4(fin,time,time_step,fout,k_tot)
call single_to_multi(fout,pop)
time = time + time_step
if (time < stop_time) goto 100
error_pop = 0.d0
do mc = 0, 1
do me = 0, 1
do mv = 0, mv_max
error_pop = error_pop + (pop(mc,me,mv) - pop_old(mc,me,mv))**2
end do
end do
end do
error_pop = SQRT(error_pop)/length_coeff
if (error_pop > max_error) then
stop_time = stop_time + 0.5d0*time_final
end if
if (error_pop > max_error) goto 100
time = 0.0d0
time = time + time_step
else
do mc = 0, 1
do me = 0, 1
do mv = 0, mv_max
pop_old(mc,me,mv) = pop(mc,me,mv)
end do
end do
end do
call multi_to_single(pop_old,fin)
call rk4(fin,time,time_step,fout,k_tot)
call single_to_multi(fout,pop)
time = time + time_step
end if
total_pop = 0.d0
mean_vib = 0.d0
do mc = 0, 1
do me = 0, 1
pop_el(mc,me) = 0.d0
end do
end do
do mv = 0, mv_max
pop_vib(mv) = 0.d0
end do
do mc = 0, 1
do me = 0, 1
do mv = 0, mv_max
total_pop = total_pop + pop(mc,me,mv)
pop_el(mc,me) = pop_el(mc,me) + pop(mc,me,mv)
end do
end do
end do
mean_vib=0.0d0
do mc = 0, 1
do me = 0, 1
do mv = 0, mv_max
mean_vib = mean_vib + mv*pop(mc,me,mv)
end do
end do
end do
i_photon = 0.d0
do mc = 0, 1
do mv = 0, mv_max
do nv = 0, mv_max
i_photon = i_photon + klu(mc,1,mv,mc,0,nv)*pop(mc,1,mv)
end do
end do
end do
current = 0.0d0
do me = 0, 1
do mv = 0, mv_max
do ne = 0, 1
do nv = 0, mv_max
current = current + klw_l(0,me,mv,1,ne,nv)*pop(0,me,mv) - klw_l(1,ne,nv,0,me,mv)*pop(1,ne,nv)
end do
end do
end do
end do
do mc = 0, 1
do me = 0, 1
do nc = 0, 1
do ne = 0, 1
c_in(mc,me,nc,ne) = 0.0d0
c_out(mc,me,nc,ne) = 0.0d0
end do
end do
end do
end do
do mv = 0, mv_max
do nv = 0, mv_max
c_in(0,0,1,0) = c_in(0,0,1,0) + klw_l(0,0,mv,1,0,nv)*pop(0,0,mv)
c_in(0,0,1,1) = c_in(0,0,1,1) + klw_l(0,0,mv,1,1,nv)*pop(0,0,mv)
c_in(0,1,1,0) = c_in(0,1,1,0) + klw_l(0,1,mv,1,0,nv)*pop(0,1,mv)
c_in(0,1,1,1) = c_in(0,1,1,1) + klw_l(0,1,mv,1,1,nv)*pop(0,1,mv)
c_out(1,0,0,0) = c_out(1,0,0,0) - klw_l(1,0,nv,0,0,mv)*pop(1,0,nv)
c_out(1,0,0,1) = c_out(1,0,0,1) - klw_l(1,0,nv,0,1,mv)*pop(1,0,nv)
c_out(1,1,0,0) = c_out(1,1,0,0) - klw_l(1,1,nv,0,0,mv)*pop(1,1,nv)
c_out(1,1,0,1) = c_out(1,1,0,1) - klw_l(1,1,nv,0,1,mv)*pop(1,1,nv)
end do
end do
!!!! Output
write (*,*) time, voltage, current*243.4d0
write (40,*) time, voltage, current*243.4d0
end do
voltage = voltage_give
if (num_voltage > 0) then voltage = voltage + (voltage_max - voltage_min)/num_voltageend if
end do
1.000000000000000E-002 40.0000000000000 1.988160357654218E-016
2.000000000000000E-002 80.0000000000000 1.364425784870382E-004
3.000000000000000E-002 120.000000000000 1.370837425057141E-004
4.000000000000000E-002 160.000000000000 2.083606552928942E-004
5.000000000000000E-002 200.000000000000 3.832251385343335E-004
6.000000000000000E-002 240.000000000000 3.846367630381155E-004
7.000000000000001E-002 280.000000000000 4.131201522155061E-004
8.000000000000000E-002 320.000000000000 5.825466767173463E-004
9.000000000000000E-002 360.000000000000 5.852940467602035E-004
0.100000000000000 400.000000000000 5.948948519956353E-004
0.110000000000000 440.000000000000 7.259287371840456E-004
0.120000000000000 480.000000000000 7.304897383455301E-004
……
38.1800000000010 2000.00000000000 74.8222446587716
9
556 Bytes, 下载次数: 7
源文件1
26.14 KB, 下载次数: 5
源文件3
5.96 KB, 下载次数: 4
源文件2
欢迎光临 Fortran Coder (http://bbs.fcode.cn/) | Powered by Discuz! X3.2 |