 implicit real*8(a-h,p-z)
 implicit integer*4(i-o)
 
 open(11,file='monomer_size.inp')
 open(12,file='potential.inp')
 open(13,file='rate_neq.inp')
 open(21,file='gpu_module.dat')
 open(22,file='cpu_module.dat')

 read(11,*) dA, dB
 read(12,*) ak_fene, r_fene, eps_vol_rep, alp_vol_rep_A, alp_vol_rep_B, alp_vol_rep_AB, eps_attr, alp_attr, delt 
 read(13,*) alambda_rep_attr, alambda_attr_neut, alambda_neut_rep

 close(11) ; close(12) ; close(13)

 pi = 4.d0*atan(1.d0)

 rA = dA/2 ; rB = dB/2

 r_fene_sq = r_fene*r_fene
 eps_fene = -0.5d0*ak_fene*r_fene_sq

 fac_vol_rep_A = 2*eps_vol_rep*alp_vol_rep_A
 fac_vol_rep_B = 2*eps_vol_rep*alp_vol_rep_B
 fac_vol_rep_AB = 2*eps_vol_rep*alp_vol_rep_AB

 rad_coro_A = 0.65d0*rA ; rad_coro_sq_A = rad_coro_A*rad_coro_A
 rad_coro_B = 0.65d0*rB ; rad_coro_sq_B = rad_coro_B*rad_coro_B
 A_kappa = 1/rA ; A_kappa_sq = A_kappa*A_kappa
 B_kappa = 1/rB ; B_kappa_sq = B_kappa*B_kappa
 fac_wall1_A =    4.0d0   *2.d0**1.5d0/rad_coro_A
 fac_wall1_B =    4.0d0   *2.d0**1.5d0/rad_coro_B
 fac_wall_ene =   4.0d0   *2.d0**1.5d0
 xi = 1.d0/(1 + 2*0.65d0*0.65d0) ; xi_2 = 2*xi ; xi_mod = xi - 0.5d0     ! kappa*rad_coro = 0.65
 zeta = sqrt(pi)*xi*erfc(0.65d0)*exp(0.65d0*0.65d0)/0.65d0 
 zeta_fac = zeta/erfc(0.65d0)
 fac_wall2 = xi_2 - 1.d0

 fac_attr = -2*eps_attr
 peak_attr = dB - 1.d0/(alp_attr*dB)

 visc = 1.d0
 fac_lang1_A = 1.d0/(6*pi*visc*rA)       ; fac_lang1_B = 1.d0/(6*pi*visc*rB)
 fac_lang2_A = sqrt(2.d0/(6*pi*visc*rA)) ; fac_lang2_B = sqrt(2.d0/(6*pi*visc*rB)) 
 pi4 = 4.d0*pi
 pi4by3 = pi4 / 3.d0
 stretch_init = 1.6d0*dB
 r_mono = rB

 del_t_sqrt = dsqrt(delt)

 write(21,*) 'real*8, constant ::  &'
 write(21,*) 'r_FPT_sq =', (2*rA)**2, ', &' 
 write(21,*) '!'
 write(21,*) 'ak_fene =', ak_fene, ', &' 
 write(21,*) 'r_fene_sq =', r_fene_sq, ', &' 
 write(21,*) 'eps_fene =', eps_fene, ', &' 
 write(21,*) '!'
 write(21,*) 'eps_vol_rep =', eps_vol_rep, ', trunc_volsq = 1.d0, &' 
 write(21,*) 'alp_vol_rep =', alp_vol_rep_A, ', &' 
 write(21,*) 'fac_vol_rep =', fac_vol_rep_A, ', &' 
 write(21,*) '!'
 write(21,*) 'eps_vol_attr =', - eps_vol_rep, ', &'
 write(21,*) 'alp_vol_attr =', alp_vol_rep_AB, ', &'
 write(21,*) 'fac_vol_attr =', -fac_vol_rep_AB, ', &'
 write(21,*) '!'
 write(21,*) 'rad_coro = ', rad_coro_A, ', &'  
 write(21,*) 'rad_coro_sq =', rad_coro_sq_A, ', &'  
 write(21,*) 'akappa =', A_kappa, ', &' 
 write(21,*) 'akappa_sq =', A_kappa_sq, ', &' 
 write(21,*) 'fac_wall1 =', fac_wall1_A, ', &'
 write(21,*) 'fac_wall_ene =', fac_wall_ene, ', &'
 write(21,*) 'xi =', xi, ', &' 
 write(21,*) 'xi_2 =', xi_2, ', &'
 write(21,*) 'xi_mod =', xi_mod, ', &'
 write(21,*) 'trunc_wall = 1.d0, &'
 write(21,*) 'zeta =', zeta, ', &'
 write(21,*) 'zeta_fac =', zeta_fac, ', &'
 write(21,*) 'fac_wall2 =', fac_wall2, ', &'
 write(21,*) '!'
 write(21,*) 'eps_attr =', eps_attr, ', &'
 write(21,*) 'alp_attr =', alp_attr, ', &'
 write(21,*) 'fac_attr =', fac_attr, ', &'
 write(21,*) 'peak_attr =', peak_attr, ', &' 
 write(21,*) '!'
 write(21,*) 'fac_lang1 =', fac_lang1_A, ', &'
 write(21,*) '!'
 write(21,*) 'alambda_rep_attr =', alambda_rep_attr, ', &'
 write(21,*) 'alambda_attr_neut =', alambda_attr_neut, ', &'
 write(21,*) 'alambda_neut_rep =', alambda_neut_rep, ', &'
 write(21,*) '!'
 write(21,*) 'delt =', delt, ', &'
 write(21,*) 'del_t_sqrt =', del_t_sqrt

 write(22,*) 'real*8, parameter ::  &'
 write(22,*) 'pi =', pi, ', &'  
 write(22,*) 'pi4 =', pi4, ', &'  
 write(22,*) 'pi4by3 =', pi4by3, ', &'  
 write(22,*) 'fac_lang2 =', fac_lang2_A, ', &'   
 write(22,*) 'stretch_init =', stretch_init, ', &'    
 write(22,*) 'r_fene =', r_fene, ', &'    
 write(22,*) 'r_mono =', r_mono 

 end
