C *******************************************************************
program canonico
implicit none
!!include 'iimc.inc'
c ** global parameters **
integer n
common /global/ n
c ** common /cofig/ Positions
doubleprecision rx(2000), ry(2000), rz(2000)
doubleprecision gdp(200)
common /cofig/ rx, ry, rz, gdp
c ** common /params/ Simulation control parameters
doubleprecision rcut, drmax, dboxmx, box
integer nstep, iratio, nblock, nrad
character disnam*20, rundat*20
common /params/ rcut, drmax, dboxmx, box,
: nstep, iratio, nblock, nrad,
: disnam, rundat
c ** common /values/ Current values of energy etc.
doubleprecision eng, pot, vir, temp, press
doubleprecision maxrad, ro
common /values/ eng, pot, vir, temp, press,
: maxrad, ro
c ** common /blkavs/ cycles average values of energy etc.
doubleprecision acp, acv, runpot, runvir
doubleprecision runrad(200)
common /blkavs/ acp, acv, runpot, runvir,
: runrad
c ** common /comran/ random numbers
integer ix1, ix2, ix3
doubleprecision rando(97)
common /comran/ rando, ix1, ix2, ix3
C *******************************************************************
C ** variable
C *******************************************************************
call config
call start
call markov
call finish
stop
end
C *******************************************************************
subroutine start
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
doubleprecision dens, dens1, vol, pots, virs, box1
doubleprecision virst, virzwz
C *******************************************************************
print *, 'iimc version 0.1'
print *, 'Montecarlo simulation of Lennard-Jones atoms'
print *, 'Periodic Boundaries in all directions'
print *, 'HECTOR DOMINGUEZ'
print *, 'Enter file to keep last running data'
read (*,'(a)') rundat
write (*,'('' file '',a)') rundat
print *,'Enter file to keep g(r) data'
read (*,'(a)') disnam
write (*,'('' file '',a)') disnam
print *, ' Enter max. displ. par. '
read *, drmax
print *, 'Enter number of blocks'
read *, nblock
print *, 'Enter number of cycles'
read *, nstep
print *, 'Enter interval for update of max. displ.'
read *, iratio
print *, 'Enter temperature'
read *, temp
print *, 'Enter potential cutoff distance'
read *, rcut
ro = box/dble(200)
nrad = int(anint(box/(2.0d0*ro)))
print *, 'No. points in distr.func. = ', nrad
print *, 'Thickness for g(r) bins = ', ro
print *, 'Number of blocks = ', nblock
print *, 'Number of cycles = ', nstep
print *, 'Ratio update interval for atoms = ', iratio
print *, 'Maximum displacement = ', drmax
print *, 'initial temperature = ', temp
print *, 'Potential cutoff distance = ', rcut
vol = box*box*box
dens = dble(n)/vol
C **** Initial configurational energy
call totener
print *, '-- Virial --=',vir
eng = 1.5d0*temp + pot/dble(n)
press = dens*temp + vir/vol
print *, 'Initial values'
print *, '# Part. = ', n
print *, 'Potential = ', pot/dble(n)
print *, 'Total eng = ', eng
print *, 'Temperature = ', temp
print *, 'Density = ', dens
print *, 'Pressure = ', press
return
end
C *******************************************************************
subroutine markov
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
integer block, step, i, j, arm, acatma, acm
doubleprecision dens, deltv, deltvi, deltvb, vol
doubleprecision ratio, dummy, ranf, nor, rno
doubleprecision vnew, vold, virnew, virold, beta
doubleprecision rxiold, ryiold, rziold, rxinew, ryinew, rzinew
doubleprecision asqpot, asqprxy, flp, flpx
C *******************************************************************
C ** Markov Chain **
write(*,'(//'' start of Markov Chain ''//)')
vol = box*box*box
dens = dble(n)/vol
beta = 1.0d0/temp
C ** Zero run accumulator **
arm = 0.0
runpot = 0.0d0
runvir = 0.0d0
asqpot = 0.0d0
asqprxy= 0.0d0
do j=1, nrad
runrad(j) = 0.0d0
enddo
write(*,'(1x,''cycle'','' Potential'','' pressure'')')
do block = 1, nblock
C ** Zero accumulators **
acatma = 0
acm = 0.0d0
acv = 0.0d0
acp = 0.0d0
C ** Particle move
do step = 1, nstep
do i = 1, n
rxiold = rx(i)
ryiold = ry(i)
rziold = rz(i)
C *** Calculate the energy of i in the old configuration ***
call enegi(rxiold,ryiold,rziold,i,vold, virold)
C *** Move i and pickup the central image ***
rxinew = rxiold + (2.0d0*ranf(dummy) - 1.0d0)*drmax
ryinew = ryiold + (2.0d0*ranf(dummy) - 1.0d0)*drmax
rzinew = rziold + (2.0d0*ranf(dummy) - 1.0d0)*drmax
rxinew = rxinew - box*anint(rxinew/box)
ryinew = ryinew - box*anint(ryinew/box)
rzinew = rzinew - box*anint(rzinew/box)
C *** Calculate the energy of i in the new configuration ***
call enegi(rxinew,ryinew,rzinew,i,vnew, virnew)
C *** Check for acceptance ***
deltv = vnew - vold
deltvi = virnew - virold
deltvb = beta*deltv
if (deltvb .lt. 75.0d0) then
if (deltvb .le. 0.0d0) then
pot = pot + deltv
vir = vir + deltvi
rx(i) = rxinew
ry(i) = ryinew
rz(i) = rzinew
acatma = acatma + 1
elseif (exp(-deltvb) .gt. ranf(dummy)) then
pot = pot + deltv
vir = vir + deltvi
rx(i) = rxinew
ry(i) = ryinew
rz(i) = rzinew
acatma = acatma + 1
endif
endif
press = dens*temp + vir/vol
C ** accumulate averages **
acm = acm + 1
acp = acp + pot/dble(n)
acv = acv + press
arm = arm + 1
runpot = runpot + pot/dble(n)
runvir = runvir + press
asqpot = asqpot + pot/dble(n)*pot/dble(n)
asqprxy= asqprxy + press*press
enddo
call radial
do j=1, nrad
runrad(j) = runrad(j)+gdp(j)
enddo
C ** Perform periodic operations **
if (mod(step,iratio) .eq. 0 ) then
C ** Adjust maximum displacement **
ratio = dble(acatma)/dble(n*iratio)
if (ratio .gt. 0.5d0) then
drmax = drmax*1.05d0
else
drmax = drmax*0.95d0
endif
acatma = 0
endif
C ** ends loop over atoms in one cycle **
enddo
nor = dble(acm)
acv = acv / nor
acp = acp / nor
eng = 1.5d0*temp + acp
C write(*,'(1x,'' ** accumulate values ** '')')
write(*,'(i5,2x,4f10.4)') block,acp,acv
enddo
C ** end loop over blocks **
C ** Normalize run accumulator **
do j=1, nrad
runrad(j) = runrad(j)/dble(nstep*nblock)
enddo
rno = dble(arm)
runpot = runpot/ rno
runvir = runvir/ rno
asqpot = (asqpot/rno) - runpot**2
asqprxy= (asqprxy/rno) - runvir**2
if (asqpot .gt. 0.0d0) flp = sqrt(asqpot)
if (asqprxy .gt. 0.0d0) flpx = sqrt(asqprxy)
eng = 1.5d0*temp + runpot
write(*,'(1x,'' ** run values ** '')')
write(*,'(1x,6x,8f10.4)') runpot,runvir
write(*,'(1x,'' ** fluctation values ** '')')
write(*,'(1x,5x,8f10.4)') flp, flpx
return
end
C *** Ends the loop over cycles ***
C *** End of Markov Chain ****
C *******************************************************************
subroutine finish
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
integer i, j
doubleprecision dens, dens1, vol
C *******************************************************************
vol = box*box*box
dens = dble(n)/vol
C **** Final total configurational energy ***
call totener
eng = 1.5d0*temp + pot/dble(n)
press = dens*temp + vir/vol
print *, 'Final values'
print *, '# Particles = ', n
print *, 'Potential = ', pot/dble(n)
print *, 'Total eng = ', eng
print *, 'Temperature = ', temp
print *, 'pressure = ', press
do i=1, n
rz(i) = rz(i) + box/2.0d0
enddo
C *** Keep data files ***
open(unit=3,file=rundat,status='unknown')
write(3,'(1x,2i14)') n
write(3,'(1x,6f14.5)') box,box,box
do i = 1, n
write(3,'(1x,3f14.5)') rx(i), ry(i), rz(i)
enddo
close(unit=3)
open(unit=3,file=disnam,status='unknown')
C write(3,'(1x,'' k'','' r'','' g(2)'')')
do i=1, nrad
write(3,'(1x,i5,2f10.5)') i,i*ro-ro/2d0,runrad(i)
enddo
close(unit=3)
return
end
C *******************************************************************
subroutine totener
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
integer i, j
doubleprecision rxi, ryi, rzi
doubleprecision rxij, ryij, rzij, rijsq, rcutsq, vcut
doubleprecision r2ij, r6ij, r12ij, fij, vij,fxi,fyi,fzi
doubleprecision fxij, fyij, fzij, boxinv
C *******************************************************************
boxinv = 1.0d0/box
rcutsq = rcut*rcut
r2ij = 1.0d0/rcutsq
r6ij = r2ij*r2ij*r2ij
r12ij = r6ij*r6ij
vcut = r12ij - r6ij
C ** Zero energies **
pot = 0.0d0
vir = 0.0d0
C ** Outer loop over i **
do i = 1, n-1
rxi = rx(i)
ryi = ry(i)
rzi = rz(i)
C ** Inner loop over j **
do j = i+1, n
rxij = rxi - rx(j)
ryij = ryi - ry(j)
rzij = rzi - rz(j)
rxij = rxij - anint(rxij*boxinv)*box
ryij = ryij - anint(ryij*boxinv)*box
rzij = rzij - anint(rzij*boxinv)*box
rijsq = rxij*rxij + ryij*ryij + rzij*rzij
if ( rijsq .lt. rcutsq ) then
r2ij = 1.0d0/rijsq
r6ij = r2ij*r2ij*r2ij
r12ij= r6ij*r6ij
vij = r12ij - r6ij
fij = (vij + r12ij)*r2ij
vij = vij - vcut
pot = pot + vij
fxij = fij*rxij
fyij = fij*ryij
fzij = fij*rzij
vir = vir + rxij*fxij + ryij*fyij + rzij*fzij
endif
enddo
enddo
C ** Incorporate correct numerical factors ***
pot = 4.0d0*pot
vir = 24.0d0*vir/3.0d0
return
end
C *******************************************************************
subroutine enegi(rxi,ryi,rzi,i,poti,viri)
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
C *******************************************************************
integer i, j
doubleprecision rxi, ryi, rzi
doubleprecision rxij, ryij, rzij, rijsq, rcutsq, vcut
doubleprecision r2ij, r6ij, r12ij, fij, vij,fxi,fyi,fzi
doubleprecision fxij, fyij, fzij, boxinv, poti,viri
C *******************************************************************
boxinv = 1.0d0/box
rcutsq = rcut*rcut
r2ij = 1.0d0/rcutsq
r6ij = r2ij*r2ij*r2ij
r12ij = r6ij*r6ij
vcut = r12ij - r6ij
C ** Zero energies **
poti = 0.0d0
viri = 0.0d0
C ** Outer loop over i **
C ** Inner loop over j **
do j = 1, n
if (j .ne. i) then
rxij = rxi - rx(j)
ryij = ryi - ry(j)
rzij = rzi - rz(j)
rxij = rxij - anint(rxij*boxinv)*box
ryij = ryij - anint(ryij*boxinv)*box
rzij = rzij - anint(rzij*boxinv)*box
rijsq = rxij*rxij + ryij*ryij + rzij*rzij
if ( rijsq .lt. rcutsq ) then
r2ij = 1.0d0/rijsq
r6ij = r2ij*r2ij*r2ij
r12ij= r6ij*r6ij
vij = r12ij - r6ij
fij = (vij + r12ij)*r2ij
vij = vij - vcut
poti = poti + vij
fxij = fij*rxij
fyij = fij*ryij
fzij = fij*rzij
viri = viri + rxij*fxij + ryij*fyij + rzij*fzij
endif
endif
enddo
C ** Incorporate correct numerical factors ***
poti = 4.0d0*poti
viri = 24.0d0*viri/3.0d0
return
end
C *******************************************************************
subroutine radial
implicit none
!!include 'iimc.inc'
C *******************************************************************
integer i, k, m, j
doubleprecision xij, yij, zij, rsq, pi, rij, boxinv
doubleprecision cte, rlower, rupper,nideal
parameter (pi=3.14159265)
C *******************************************************************
boxinv = 1.0d0/box
maxrad = box/2.0d0
ro = maxrad/dble(nrad)
do j=1, nrad
gdp(j) = 0.0d0
enddo
do i=1, n-1
do j = i+1, n
xij = rx(i)-rx(j)
yij = ry(i)-ry(j)
zij = rz(i)-rz(j)
xij = xij - box*anint(xij*boxinv)
yij = yij - box*anint(yij*boxinv)
zij = zij - box*anint(zij*boxinv)
rsq = xij*xij + yij*yij + zij*zij
rij = sqrt(rsq)
k = int(rij/ro)+1
if (k .le. nrad) then
gdp(k) = gdp(k)+2
endif
enddo
enddo
cte = 4.0d0*pi*(dble(n)/(box*box*box))/3.0d0
do m=1, nrad
rlower = dble(m-1)*ro
rupper = rlower+ro
nideal = cte*(rupper**3 - rlower**3)
gdp(m) = gdp(m)/(dble(n)*nideal)
enddo
return
end
C *******************************************************************
subroutine config
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
integer i
character name*20
C *******************************************************************
print *,'enter file'
read (*,'(a)') name
write(*,'('' file '',a)') name
C ******* read initial configuration *******
open (unit=3, file= name, status= 'old')
read(3,*) n
read(3,*) box, box, box
do i=1, n
read(3,*) rx(i), ry(i), rz(i)
enddo
close (unit=3)
do i=1, n
rz(i) = rz(i) - box/2.0d0
enddo
C ** random number ***
call ranset
return
end
C *******************************************************************
subroutine ranset
implicit none
!!include 'iimc.inc'
C ********************************************************************
C ***** local variables *****
integer m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3
integer j, idum, seed
doubleprecision rm1, rm2
parameter (m1=259200, ia1=7141, ic1=54773)
parameter (m2=134456, ia2=8121, ic2=28411)
parameter (m3=243000, ia3=4561, ic3=51349)
parameter (rm1=1.0/m1, rm2=1.0/m2)
C ********************************************************************
print *, 'Enter random number integer'
read *, seed
idum = 2.0*seed + 1.0
ix1 = mod(ic1+idum,m1)
ix1 = mod(ia1*ix1+ic1,m1)
ix2 = mod(ix1,m2)
ix1 = mod(ia1*ix1+ic1,m1)
ix3 = mod(ix1,m3)
do j=1, 97
ix1 = mod(ia1*ix1+ic1,m1)
ix1 = mod(ia2*ix2+ic2,m2)
rando(j) = (float(ix1) + float(ix2)*rm2)*rm1
enddo
return
end
C *******************************************************************
doubleprecision function ranf (dummy)
implicit none
!!include 'iimc.inc'
C ********************************************************************
C ***** local variables *****
integer m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3, j
doubleprecision rm1, rm2, dummy
parameter (m1=259200, ia1=7141, ic1=54773)
parameter (m2=134456, ia2=8121, ic2=28411)
parameter (m3=243000, ia3=4561, ic3=51349)
parameter (rm1=1.0/m1, rm2=1.0/m2)
C ********************************************************************
10 ix1 = mod(ia1*ix1+ic1,m1)
ix2 = mod(ia2*ix2+ic2,m2)
ix3 = mod(ia3*ix3+ic3,m3)
20 j = 1 + (97*ix3)/m3
if ((j .gt. 97) .or. (j .lt. 1)) then
print *, 'j error', j
goto 20
endif
ranf = rando(j)
rando(j) = (float(ix1) + float(ix2)*rm2)*rm1
if ((ranf .lt. 0.0) .or. (ranf .gt. 1)) then
print *, 'error in random number'
goto 10
endif
return
end
end return endif goto 10 print *, 'error in random number' if ((ranf .lt. 0.0) .or. (ranf .gt. 1)) then rando(j) = (float(ix1) + float(ix2)*rm2)*rm1 ranf = rando(j) endif goto 20 print *, 'j error', j if ((j .gt. 97) .or. (j .lt. 1)) then20 j = 1 + (97*ix3)/m3 ix3 = mod(ia3*ix3+ic3,m3) ix2 = mod(ia2*ix2+ic2,m2)10 ix1 = mod(ia1*ix1+ic1,m1)C ******************************************************************** parameter (rm1=1.0/m1, rm2=1.0/m2) parameter (m3=243000, ia3=4561, ic3=51349) parameter (m2=134456, ia2=8121, ic2=28411) parameter (m1=259200, ia1=7141, ic1=54773) doubleprecision rm1, rm2, dummy integer m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3, jC ***** local variables ***** C ******************************************************************** !!include 'iimc.inc' implicit none doubleprecision function ranf (dummy)C ******************************************************************* end return enddo rando(j) = (float(ix1) + float(ix2)*rm2)*rm1 ix1 = mod(ia2*ix2+ic2,m2) ix1 = mod(ia1*ix1+ic1,m1) do j=1, 97 ix3 = mod(ix1,m3) ix1 = mod(ia1*ix1+ic1,m1) ix2 = mod(ix1,m2) ix1 = mod(ia1*ix1+ic1,m1) ix1 = mod(ic1+idum,m1) idum = 2.0*seed + 1.0 read *, seed print *, 'Enter random number integer'C ******************************************************************** parameter (rm1=1.0/m1, rm2=1.0/m2) parameter (m3=243000, ia3=4561, ic3=51349) parameter (m2=134456, ia2=8121, ic2=28411) parameter (m1=259200, ia1=7141, ic1=54773) doubleprecision rm1, rm2 integer j, idum, seed integer m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3C ***** local variables ***** C ******************************************************************** !!include 'iimc.inc' implicit none subroutine ransetC ******************************************************************* end return call ransetC ** random number *** enddo rz(i) = rz(i) - box/2.0d0 do i=1, n close (unit=3) enddo read(3,*) rx(i), ry(i), rz(i) do i=1, n read(3,*) box, box, box read(3,*) n open (unit=3, file= name, status= 'old')C ******* read initial configuration ******* write(*,'('' file '',a)') name read (*,'(a)') name print *,'enter file'C ******************************************************************* character name*20 integer iC ** Local variables **C ******************************************************************* !!include 'iimc.inc' implicit none subroutine configC ******************************************************************* end return enddo gdp(m) = gdp(m)/(dble(n)*nideal) nideal = cte*(rupper**3 - rlower**3) rupper = rlower+ro rlower = dble(m-1)*ro do m=1, nrad cte = 4.0d0*pi*(dble(n)/(box*box*box))/3.0d0 enddo enddo endif gdp(k) = gdp(k)+2 if (k .le. nrad) then k = int(rij/ro)+1 rij = sqrt(rsq) rsq = xij*xij + yij*yij + zij*zij zij = zij - box*anint(zij*boxinv) yij = yij - box*anint(yij*boxinv) xij = xij - box*anint(xij*boxinv) zij = rz(i)-rz(j) yij = ry(i)-ry(j) xij = rx(i)-rx(j) do j = i+1, n do i=1, n-1 enddo gdp(j) = 0.0d0 do j=1, nrad ro = maxrad/dble(nrad) maxrad = box/2.0d0 boxinv = 1.0d0/boxC ******************************************************************* parameter (pi=3.14159265) doubleprecision cte, rlower, rupper,nideal doubleprecision xij, yij, zij, rsq, pi, rij, boxinv integer i, k, m, jC ******************************************************************* !!include 'iimc.inc' implicit none subroutine radialC ******************************************************************* end return viri = 24.0d0*viri/3.0d0 poti = 4.0d0*potiC ** Incorporate correct numerical factors *** enddo endif endif viri = viri + rxij*fxij + ryij*fyij + rzij*fzij fzij = fij*rzij fyij = fij*ryij fxij = fij*rxij poti = poti + vij vij = vij - vcut fij = (vij + r12ij)*r2ij vij = r12ij - r6ij r12ij= r6ij*r6ij r6ij = r2ij*r2ij*r2ij r2ij = 1.0d0/rijsq if ( rijsq .lt. rcutsq ) then rijsq = rxij*rxij + ryij*ryij + rzij*rzij rzij = rzij - anint(rzij*boxinv)*box ryij = ryij - anint(ryij*boxinv)*box rxij = rxij - anint(rxij*boxinv)*box rzij = rzi - rz(j) ryij = ryi - ry(j) rxij = rxi - rx(j) if (j .ne. i) then do j = 1, nC ** Inner loop over j **C ** Outer loop over i ** viri = 0.0d0 poti = 0.0d0C ** Zero energies ** vcut = r12ij - r6ij r12ij = r6ij*r6ij r6ij = r2ij*r2ij*r2ij r2ij = 1.0d0/rcutsq rcutsq = rcut*rcut boxinv = 1.0d0/boxC ******************************************************************* doubleprecision fxij, fyij, fzij, boxinv, poti,viri doubleprecision r2ij, r6ij, r12ij, fij, vij,fxi,fyi,fzi doubleprecision rxij, ryij, rzij, rijsq, rcutsq, vcut doubleprecision rxi, ryi, rzi integer i, jC *******************************************************************C ** Local variables **C ******************************************************************* !!include 'iimc.inc' implicit none subroutine enegi(rxi,ryi,rzi,i,poti,viri)C ******************************************************************* end return vir = 24.0d0*vir/3.0d0 pot = 4.0d0*potC ** Incorporate correct numerical factors *** enddo enddo endif vir = vir + rxij*fxij + ryij*fyij + rzij*fzij fzij = fij*rzij fyij = fij*ryij fxij = fij*rxij pot = pot + vij vij = vij - vcut fij = (vij + r12ij)*r2ij vij = r12ij - r6ij r12ij= r6ij*r6ij r6ij = r2ij*r2ij*r2ij r2ij = 1.0d0/rijsq if ( rijsq .lt. rcutsq ) then rijsq = rxij*rxij + ryij*ryij + rzij*rzij rzij = rzij - anint(rzij*boxinv)*box ryij = ryij - anint(ryij*boxinv)*box rxij = rxij - anint(rxij*boxinv)*box rzij = rzi - rz(j) ryij = ryi - ry(j) rxij = rxi - rx(j) do j = i+1, nC ** Inner loop over j ** rzi = rz(i) ryi = ry(i) rxi = rx(i) do i = 1, n-1 C ** Outer loop over i ** vir = 0.0d0 pot = 0.0d0C ** Zero energies ** vcut = r12ij - r6ij r12ij = r6ij*r6ij r6ij = r2ij*r2ij*r2ij r2ij = 1.0d0/rcutsq rcutsq = rcut*rcut boxinv = 1.0d0/boxC ******************************************************************* doubleprecision fxij, fyij, fzij, boxinv doubleprecision r2ij, r6ij, r12ij, fij, vij,fxi,fyi,fzi doubleprecision rxij, ryij, rzij, rijsq, rcutsq, vcut doubleprecision rxi, ryi, rzi integer i, jC ** Local variables **C ******************************************************************* !!include 'iimc.inc' implicit none subroutine totenerC ******************************************************************* end return close(unit=3) enddo write(3,'(1x,i5,2f10.5)') i,i*ro-ro/2d0,runrad(i) do i=1, nradC write(3,'(1x,'' k'','' r'','' g(2)'')') open(unit=3,file=disnam,status='unknown') close(unit=3) enddo write(3,'(1x,3f14.5)') rx(i), ry(i), rz(i) do i = 1, n write(3,'(1x,6f14.5)') box,box,box write(3,'(1x,2i14)') n open(unit=3,file=rundat,status='unknown')C *** Keep data files *** enddo rz(i) = rz(i) + box/2.0d0 do i=1, n print *, 'pressure = ', press print *, 'Temperature = ', temp print *, 'Total eng = ', eng print *, 'Potential = ', pot/dble(n) print *, '# Particles = ', n print *, 'Final values' press = dens*temp + vir/vol eng = 1.5d0*temp + pot/dble(n) call totener C **** Final total configurational energy *** dens = dble(n)/vol vol = box*box*boxC ******************************************************************* doubleprecision dens, dens1, vol integer i, jC ** Local variables **C ******************************************************************* !!include 'iimc.inc' implicit none subroutine finishC ******************************************************************* C *** End of Markov Chain ****C *** Ends the loop over cycles *** end return write(*,'(1x,5x,8f10.4)') flp, flpx write(*,'(1x,'' ** fluctation values ** '')') write(*,'(1x,6x,8f10.4)') runpot,runvir write(*,'(1x,'' ** run values ** '')') eng = 1.5d0*temp + runpot if (asqprxy .gt. 0.0d0) flpx = sqrt(asqprxy) if (asqpot .gt. 0.0d0) flp = sqrt(asqpot) asqprxy= (asqprxy/rno) - runvir**2 asqpot = (asqpot/rno) - runpot**2 runvir = runvir/ rno runpot = runpot/ rno rno = dble(arm) enddo runrad(j) = runrad(j)/dble(nstep*nblock) do j=1, nradC ** Normalize run accumulator **C ** end loop over blocks ** enddo write(*,'(i5,2x,4f10.4)') block,acp,acvC write(*,'(1x,'' ** accumulate values ** '')') eng = 1.5d0*temp + acp acp = acp / nor acv = acv / nor nor = dble(acm) enddoC ** ends loop over atoms in one cycle ** endif acatma = 0 endif drmax = drmax*0.95d0 else drmax = drmax*1.05d0 if (ratio .gt. 0.5d0) then ratio = dble(acatma)/dble(n*iratio)C ** Adjust maximum displacement ** if (mod(step,iratio) .eq. 0 ) thenC ** Perform periodic operations ** enddo runrad(j) = runrad(j)+gdp(j) do j=1, nrad call radial enddo asqprxy= asqprxy + press*press asqpot = asqpot + pot/dble(n)*pot/dble(n) runvir = runvir + press runpot = runpot + pot/dble(n) arm = arm + 1 acv = acv + press acp = acp + pot/dble(n) acm = acm + 1C ** accumulate averages ** press = dens*temp + vir/vol endif endif acatma = acatma + 1 rz(i) = rzinew ry(i) = ryinew rx(i) = rxinew vir = vir + deltvi pot = pot + deltv elseif (exp(-deltvb) .gt. ranf(dummy)) then acatma = acatma + 1 rz(i) = rzinew ry(i) = ryinew rx(i) = rxinew vir = vir + deltvi pot = pot + deltv if (deltvb .le. 0.0d0) then if (deltvb .lt. 75.0d0) then deltvb = beta*deltv deltvi = virnew - virold deltv = vnew - vold C *** Check for acceptance *** call enegi(rxinew,ryinew,rzinew,i,vnew, virnew)C *** Calculate the energy of i in the new configuration *** rzinew = rzinew - box*anint(rzinew/box) ryinew = ryinew - box*anint(ryinew/box) rxinew = rxinew - box*anint(rxinew/box) rzinew = rziold + (2.0d0*ranf(dummy) - 1.0d0)*drmax ryinew = ryiold + (2.0d0*ranf(dummy) - 1.0d0)*drmax rxinew = rxiold + (2.0d0*ranf(dummy) - 1.0d0)*drmaxC *** Move i and pickup the central image *** call enegi(rxiold,ryiold,rziold,i,vold, virold) C *** Calculate the energy of i in the old configuration *** rziold = rz(i) ryiold = ry(i) rxiold = rx(i) do i = 1, n do step = 1, nstepC ** Particle move acp = 0.0d0 acv = 0.0d0 acm = 0.0d0 acatma = 0C ** Zero accumulators ** do block = 1, nblock write(*,'(1x,''cycle'','' Potential'','' pressure'')') enddo runrad(j) = 0.0d0 do j=1, nrad asqprxy= 0.0d0 asqpot = 0.0d0 runvir = 0.0d0 runpot = 0.0d0 arm = 0.0C ** Zero run accumulator ** beta = 1.0d0/temp dens = dble(n)/vol vol = box*box*box write(*,'(//'' start of Markov Chain ''//)')C ** Markov Chain **C ******************************************************************* doubleprecision asqpot, asqprxy, flp, flpx doubleprecision rxiold, ryiold, rziold, rxinew, ryinew, rzinew doubleprecision vnew, vold, virnew, virold, beta doubleprecision ratio, dummy, ranf, nor, rno doubleprecision dens, deltv, deltvi, deltvb, vol integer block, step, i, j, arm, acatma, acmC ** Local variables **C ******************************************************************* !!include 'iimc.inc' implicit none subroutine markovC ******************************************************************* end return print *, 'Pressure = ', press print *, 'Density = ', dens print *, 'Temperature = ', temp print *, 'Total eng = ', eng print *, 'Potential = ', pot/dble(n) print *, '# Part. = ', n print *, 'Initial values' press = dens*temp + vir/vol eng = 1.5d0*temp + pot/dble(n) print *, '-- Virial --=',vir call totener C **** Initial configurational energy dens = dble(n)/vol vol = box*box*box print *, 'Potential cutoff distance = ', rcut print *, 'initial temperature = ', temp print *, 'Maximum displacement = ', drmax print *, 'Ratio update interval for atoms = ', iratio print *, 'Number of cycles = ', nstep print *, 'Number of blocks = ', nblock print *, 'Thickness for g(r) bins = ', ro print *, 'No. points in distr.func. = ', nrad nrad = int(anint(box/(2.0d0*ro))) ro = box/dble(200) read *, rcut print *, 'Enter potential cutoff distance' read *, temp print *, 'Enter temperature' read *, iratio print *, 'Enter interval for update of max. displ.' read *, nstep print *, 'Enter number of cycles' read *, nblock print *, 'Enter number of blocks' read *, drmax print *, ' Enter max. displ. par. ' write (*,'('' file '',a)') disnam read (*,'(a)') disnam print *,'Enter file to keep g(r) data' write (*,'('' file '',a)') rundat read (*,'(a)') rundat print *, 'Enter file to keep last running data' print *, 'HECTOR DOMINGUEZ' print *, 'Periodic Boundaries in all directions' print *, 'Montecarlo simulation of Lennard-Jones atoms' print *, 'iimc version 0.1' C ******************************************************************* doubleprecision virst, virzwz doubleprecision dens, dens1, vol, pots, virs, box1C ** Local variables **C ******************************************************************* !!include 'iimc.inc' implicit none subroutine startC ******************************************************************* end stop call finish call markov call start call configC *******************************************************************C ** variable C ******************************************************************* !!include 'iimc.inc' implicit none program canonicoC *******************************************************************
C *******************************************************************
program canonico
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** variable
C *******************************************************************
call config
call start
call markov
call finish
stop
end
C *******************************************************************
subroutine start
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
doubleprecision dens, dens1, vol, pots, virs, box1
doubleprecision virst, virzwz
C *******************************************************************
print *, 'iimc version 0.1'
print *, 'Montecarlo simulation of Lennard-Jones atoms'
print *, 'Periodic Boundaries in all directions'
print *, 'HECTOR DOMINGUEZ'
print *, 'Enter file to keep last running data'
read (*,'(a)') rundat
write (*,'('' file '',a)') rundat
print *,'Enter file to keep g(r) data'
read (*,'(a)') disnam
write (*,'('' file '',a)') disnam
print *, ' Enter max. displ. par. '
read *, drmax
print *, 'Enter number of blocks'
read *, nblock
print *, 'Enter number of cycles'
read *, nstep
print *, 'Enter interval for update of max. displ.'
read *, iratio
print *, 'Enter temperature'
read *, temp
print *, 'Enter potential cutoff distance'
read *, rcut
ro = box/dble(200)
nrad = int(anint(box/(2.0d0*ro)))
print *, 'No. points in distr.func. = ', nrad
print *, 'Thickness for g(r) bins = ', ro
print *, 'Number of blocks = ', nblock
print *, 'Number of cycles = ', nstep
print *, 'Ratio update interval for atoms = ', iratio
print *, 'Maximum displacement = ', drmax
print *, 'initial temperature = ', temp
print *, 'Potential cutoff distance = ', rcut
vol = box*box*box
dens = dble(n)/vol
C **** Initial configurational energy
call totener
print *, '-- Virial --=',vir
eng = 1.5d0*temp + pot/dble(n)
press = dens*temp + vir/vol
print *, 'Initial values'
print *, '# Part. = ', n
print *, 'Potential = ', pot/dble(n)
print *, 'Total eng = ', eng
print *, 'Temperature = ', temp
print *, 'Density = ', dens
print *, 'Pressure = ', press
return
end
C *******************************************************************
subroutine markov
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
integer block, step, i, j, arm, acatma, acm
doubleprecision dens, deltv, deltvi, deltvb, vol
doubleprecision ratio, dummy, ranf, nor, rno
doubleprecision vnew, vold, virnew, virold, beta
doubleprecision rxiold, ryiold, rziold, rxinew, ryinew, rzinew
doubleprecision asqpot, asqprxy, flp, flpx
C *******************************************************************
C ** Markov Chain **
write(*,'(//'' start of Markov Chain ''//)')
vol = box*box*box
dens = dble(n)/vol
beta = 1.0d0/temp
C ** Zero run accumulator **
arm = 0.0
runpot = 0.0d0
runvir = 0.0d0
asqpot = 0.0d0
asqprxy= 0.0d0
do j=1, nrad
runrad(j) = 0.0d0
enddo
write(*,'(1x,''cycle'','' Potential'','' pressure'')')
do block = 1, nblock
C ** Zero accumulators **
acatma = 0
acm = 0.0d0
acv = 0.0d0
acp = 0.0d0
C ** Particle move
do step = 1, nstep
do i = 1, n
rxiold = rx(i)
ryiold = ry(i)
rziold = rz(i)
C *** Calculate the energy of i in the old configuration ***
call enegi(rxiold,ryiold,rziold,i,vold, virold)
C *** Move i and pickup the central image ***
rxinew = rxiold + (2.0d0*ranf(dummy) - 1.0d0)*drmax
ryinew = ryiold + (2.0d0*ranf(dummy) - 1.0d0)*drmax
rzinew = rziold + (2.0d0*ranf(dummy) - 1.0d0)*drmax
rxinew = rxinew - box*anint(rxinew/box)
ryinew = ryinew - box*anint(ryinew/box)
rzinew = rzinew - box*anint(rzinew/box)
C *** Calculate the energy of i in the new configuration ***
call enegi(rxinew,ryinew,rzinew,i,vnew, virnew)
C *** Check for acceptance ***
deltv = vnew - vold
deltvi = virnew - virold
deltvb = beta*deltv
if (deltvb .lt. 75.0d0) then
if (deltvb .le. 0.0d0) then
pot = pot + deltv
vir = vir + deltvi
rx(i) = rxinew
ry(i) = ryinew
rz(i) = rzinew
acatma = acatma + 1
elseif (exp(-deltvb) .gt. ranf(dummy)) then
pot = pot + deltv
vir = vir + deltvi
rx(i) = rxinew
ry(i) = ryinew
rz(i) = rzinew
acatma = acatma + 1
endif
endif
press = dens*temp + vir/vol
C ** accumulate averages **
acm = acm + 1
acp = acp + pot/dble(n)
acv = acv + press
arm = arm + 1
runpot = runpot + pot/dble(n)
runvir = runvir + press
asqpot = asqpot + pot/dble(n)*pot/dble(n)
asqprxy= asqprxy + press*press
enddo
call radial
do j=1, nrad
runrad(j) = runrad(j)+gdp(j)
enddo
C ** Perform periodic operations **
if (mod(step,iratio) .eq. 0 ) then
C ** Adjust maximum displacement **
ratio = dble(acatma)/dble(n*iratio)
if (ratio .gt. 0.5d0) then
drmax = drmax*1.05d0
else
drmax = drmax*0.95d0
endif
acatma = 0
endif
C ** ends loop over atoms in one cycle **
enddo
nor = dble(acm)
acv = acv / nor
acp = acp / nor
eng = 1.5d0*temp + acp
C write(*,'(1x,'' ** accumulate values ** '')')
write(*,'(i5,2x,4f10.4)') block,acp,acv
enddo
C ** end loop over blocks **
C ** Normalize run accumulator **
do j=1, nrad
runrad(j) = runrad(j)/dble(nstep*nblock)
enddo
rno = dble(arm)
runpot = runpot/ rno
runvir = runvir/ rno
asqpot = (asqpot/rno) - runpot**2
asqprxy= (asqprxy/rno) - runvir**2
if (asqpot .gt. 0.0d0) flp = sqrt(asqpot)
if (asqprxy .gt. 0.0d0) flpx = sqrt(asqprxy)
eng = 1.5d0*temp + runpot
write(*,'(1x,'' ** run values ** '')')
write(*,'(1x,6x,8f10.4)') runpot,runvir
write(*,'(1x,'' ** fluctation values ** '')')
write(*,'(1x,5x,8f10.4)') flp, flpx
return
end
C *** Ends the loop over cycles ***
C *** End of Markov Chain ****
C *******************************************************************
subroutine finish
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
integer i, j
doubleprecision dens, dens1, vol
C *******************************************************************
vol = box*box*box
dens = dble(n)/vol
C **** Final total configurational energy ***
call totener
eng = 1.5d0*temp + pot/dble(n)
press = dens*temp + vir/vol
print *, 'Final values'
print *, '# Particles = ', n
print *, 'Potential = ', pot/dble(n)
print *, 'Total eng = ', eng
print *, 'Temperature = ', temp
print *, 'pressure = ', press
do i=1, n
rz(i) = rz(i) + box/2.0d0
enddo
C *** Keep data files ***
open(unit=3,file=rundat,status='unknown')
write(3,'(1x,2i14)') n
write(3,'(1x,6f14.5)') box,box,box
do i = 1, n
write(3,'(1x,3f14.5)') rx(i), ry(i), rz(i)
enddo
close(unit=3)
open(unit=3,file=disnam,status='unknown')
C write(3,'(1x,'' k'','' r'','' g(2)'')')
do i=1, nrad
write(3,'(1x,i5,2f10.5)') i,i*ro-ro/2d0,runrad(i)
enddo
close(unit=3)
return
end
C *******************************************************************
subroutine totener
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
integer i, j
doubleprecision rxi, ryi, rzi
doubleprecision rxij, ryij, rzij, rijsq, rcutsq, vcut
doubleprecision r2ij, r6ij, r12ij, fij, vij,fxi,fyi,fzi
doubleprecision fxij, fyij, fzij, boxinv
C *******************************************************************
boxinv = 1.0d0/box
rcutsq = rcut*rcut
r2ij = 1.0d0/rcutsq
r6ij = r2ij*r2ij*r2ij
r12ij = r6ij*r6ij
vcut = r12ij - r6ij
C ** Zero energies **
pot = 0.0d0
vir = 0.0d0
C ** Outer loop over i **
do i = 1, n-1
rxi = rx(i)
ryi = ry(i)
rzi = rz(i)
C ** Inner loop over j **
do j = i+1, n
rxij = rxi - rx(j)
ryij = ryi - ry(j)
rzij = rzi - rz(j)
rxij = rxij - anint(rxij*boxinv)*box
ryij = ryij - anint(ryij*boxinv)*box
rzij = rzij - anint(rzij*boxinv)*box
rijsq = rxij*rxij + ryij*ryij + rzij*rzij
if ( rijsq .lt. rcutsq ) then
r2ij = 1.0d0/rijsq
r6ij = r2ij*r2ij*r2ij
r12ij= r6ij*r6ij
vij = r12ij - r6ij
fij = (vij + r12ij)*r2ij
vij = vij - vcut
pot = pot + vij
fxij = fij*rxij
fyij = fij*ryij
fzij = fij*rzij
vir = vir + rxij*fxij + ryij*fyij + rzij*fzij
endif
enddo
enddo
C ** Incorporate correct numerical factors ***
pot = 4.0d0*pot
vir = 24.0d0*vir/3.0d0
return
end
C *******************************************************************
subroutine enegi(rxi,ryi,rzi,i,poti,viri)
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
C *******************************************************************
integer i, j
doubleprecision rxi, ryi, rzi
doubleprecision rxij, ryij, rzij, rijsq, rcutsq, vcut
doubleprecision r2ij, r6ij, r12ij, fij, vij,fxi,fyi,fzi
doubleprecision fxij, fyij, fzij, boxinv, poti,viri
C *******************************************************************
boxinv = 1.0d0/box
rcutsq = rcut*rcut
r2ij = 1.0d0/rcutsq
r6ij = r2ij*r2ij*r2ij
r12ij = r6ij*r6ij
vcut = r12ij - r6ij
C ** Zero energies **
poti = 0.0d0
viri = 0.0d0
C ** Outer loop over i **
C ** Inner loop over j **
do j = 1, n
if (j .ne. i) then
rxij = rxi - rx(j)
ryij = ryi - ry(j)
rzij = rzi - rz(j)
rxij = rxij - anint(rxij*boxinv)*box
ryij = ryij - anint(ryij*boxinv)*box
rzij = rzij - anint(rzij*boxinv)*box
rijsq = rxij*rxij + ryij*ryij + rzij*rzij
if ( rijsq .lt. rcutsq ) then
r2ij = 1.0d0/rijsq
r6ij = r2ij*r2ij*r2ij
r12ij= r6ij*r6ij
vij = r12ij - r6ij
fij = (vij + r12ij)*r2ij
vij = vij - vcut
poti = poti + vij
fxij = fij*rxij
fyij = fij*ryij
fzij = fij*rzij
viri = viri + rxij*fxij + ryij*fyij + rzij*fzij
endif
endif
enddo
C ** Incorporate correct numerical factors ***
poti = 4.0d0*poti
viri = 24.0d0*viri/3.0d0
return
end
C *******************************************************************
subroutine radial
implicit none
!!include 'iimc.inc'
C *******************************************************************
integer i, k, m, j
doubleprecision xij, yij, zij, rsq, pi, rij, boxinv
doubleprecision cte, rlower, rupper,nideal
parameter (pi=3.14159265)
C *******************************************************************
boxinv = 1.0d0/box
maxrad = box/2.0d0
ro = maxrad/dble(nrad)
do j=1, nrad
gdp(j) = 0.0d0
enddo
do i=1, n-1
do j = i+1, n
xij = rx(i)-rx(j)
yij = ry(i)-ry(j)
zij = rz(i)-rz(j)
xij = xij - box*anint(xij*boxinv)
yij = yij - box*anint(yij*boxinv)
zij = zij - box*anint(zij*boxinv)
rsq = xij*xij + yij*yij + zij*zij
rij = sqrt(rsq)
k = int(rij/ro)+1
if (k .le. nrad) then
gdp(k) = gdp(k)+2
endif
enddo
enddo
cte = 4.0d0*pi*(dble(n)/(box*box*box))/3.0d0
do m=1, nrad
rlower = dble(m-1)*ro
rupper = rlower+ro
nideal = cte*(rupper**3 - rlower**3)
gdp(m) = gdp(m)/(dble(n)*nideal)
enddo
return
end
C *******************************************************************
subroutine config
implicit none
!!include 'iimc.inc'
C *******************************************************************
C ** Local variables **
integer i
character name*20
C *******************************************************************
print *,'enter file'
read (*,'(a)') name
write(*,'('' file '',a)') name
C ******* read initial configuration *******
open (unit=3, file= name, status= 'old')
read(3,*) n
read(3,*) box, box, box
do i=1, n
read(3,*) rx(i), ry(i), rz(i)
enddo
close (unit=3)
do i=1, n
rz(i) = rz(i) - box/2.0d0
enddo
C ** random number ***
call ranset
return
end
C *******************************************************************
subroutine ranset
implicit none
!!include 'iimc.inc'
C ********************************************************************
C ***** local variables *****
integer m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3
integer j, idum, seed
doubleprecision rm1, rm2
parameter (m1=259200, ia1=7141, ic1=54773)
parameter (m2=134456, ia2=8121, ic2=28411)
parameter (m3=243000, ia3=4561, ic3=51349)
parameter (rm1=1.0/m1, rm2=1.0/m2)
C ********************************************************************
print *, 'Enter random number integer'
read *, seed
idum = 2.0*seed + 1.0
ix1 = mod(ic1+idum,m1)
ix1 = mod(ia1*ix1+ic1,m1)
ix2 = mod(ix1,m2)
ix1 = mod(ia1*ix1+ic1,m1)
ix3 = mod(ix1,m3)
do j=1, 97
ix1 = mod(ia1*ix1+ic1,m1)
ix1 = mod(ia2*ix2+ic2,m2)
rando(j) = (float(ix1) + float(ix2)*rm2)*rm1
enddo
return
end
C *******************************************************************
doubleprecision function ranf (dummy)
implicit none
!!include 'iimc.inc'
C ********************************************************************
C ***** local variables *****
integer m1, m2, m3, ia1, ia2, ia3, ic1, ic2, ic3, j
doubleprecision rm1, rm2, dummy
parameter (m1=259200, ia1=7141, ic1=54773)
parameter (m2=134456, ia2=8121, ic2=28411)
parameter (m3=243000, ia3=4561, ic3=51349)
parameter (rm1=1.0/m1, rm2=1.0/m2)
C ********************************************************************
10 ix1 = mod(ia1*ix1+ic1,m1)
ix2 = mod(ia2*ix2+ic2,m2)
ix3 = mod(ia3*ix3+ic3,m3)
20 j = 1 + (97*ix3)/m3
if ((j .gt. 97) .or. (j .lt. 1)) then
print *, 'j error', j
goto 20
endif
ranf = rando(j)
rando(j) = (float(ix1) + float(ix2)*rm2)*rm1
if ((ranf .lt. 0.0) .or. (ranf .gt. 1)) then
print *, 'error in random number'
goto 10
endif
return
end
program map
integer*4 i,imax,j
! imax gives maximum number of points per orbit
parameter(imax=1000)
real*8 x,y,xnew,ynew,alpha,beta,r
real*8 xdata(imax),ydata(imax)
open(11, file='map_data')
! get alpha
write(*,*) 'give value of alpha,beta'
read(*,*) alpha,beta
write(*,*) alpha,beta
! loop over 101 starting points to build up graph
do j=0,100
! initialize i and x,y
i=1
x = 0
y = j * 0.01
xdata(1)=x
ydata(1)=y
! start follow trajectory
do while (i .lt. imax )
! calculate new point
xnew = y + 1 - alpha * x * x
ynew = beta * x
! calculate distance square from origin
r = xnew*xnew+ynew*ynew
if ( r .lt. 100 ) then
! if not running away store point and update
x = xnew
y = ynew
i=i+1
xdata(i)=x
ydata(i)=y
else
! this is a diverging orbit, terminate
i=imax
xdata(imax)=10000
end if
end do
! print if not diverging trajectory
if ( xdata(imax) .lt. 100) then
write(11,1000) (xdata(i),ydata(i),i=100,imax)
end if
1000 format(2f10.5)
! xmgr did not like this one
! write(11,*)
end stop 'map' end do ! write(11,*) ! xmgr did not like this one 1000 format(2f10.5) end if write(11,1000) (xdata(i),ydata(i),i=1,imax) if ( xdata(imax) .lt. 100) then ! print if not diverging trajectory end do end if xdata(imax)=10000 i=imax ! this is a diverging orbit, terminate else ydata(i)=y xdata(i)=x i=i+1 y = ynew x = xnew ! if not running away store point and update if ( r .lt. 100 ) then r = xnew*xnew+ynew*ynew ! calculate distance square from origin ynew = x * sa + c * ca xnew = x * ca - c * sa c=y - x*x ! calculate new point do while (i .lt. imax ) ! start follow trajectory ydata(1)=y xdata(1)=x y = j * 0.01 x = 0 i=1 ! initialize i and x,y do j=0,100 ! loop over 101 starting points to build up graph sa=sin(alpha) ca=cos(alpha) ! calculate cos and sin read(*,*) alpha write(*,*) 'give value of alpha' ! get alpha open(11, file='map_data') real*8 xdata(imax),ydata(imax) real*8 x,y,xnew,ynew,alpha,ca,sa,c,r parameter(imax=1000) ! imax gives maximum number of points per orbit integer*4 i,imax,j program map
end do
stop 'map'
end
end return#endif $ nfftdim1,nfftdim2,fftable,nfftable,ffwork,nffwork) call pubz3d(isign,nfft1,nfft2,nfft3,array,#else $ ffwork,isys) $ nfftdim1,nfftdim2,array,nfftdim1,nfftdim2,fftable, call CCFFT3D(isign,nfft1,nfft2,nfft3,scale,array, isys(4)=0 isys(3)=0 isys(2)=0 isys(1)=3 scale = 1.d0#ifdef CRAY#endif $ nfftdim1,nfftdim2,fftable) call ZFFT3D(isign,nfft1,nfft2,nfft3,array,#ifdef SGIFFT isign = -1 double precision scale integer isign,inc1,inc2,inc3 integer nfftable,nffwork,isys(4) integer nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3 double precision array(*),fftable(*),ffwork(*) implicit none $ nfftable,nffwork) $ nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3, subroutine fft_back(array,fftable,ffwork,c----------------------------------------------------------- end return#endif $ nfftdim1,nfftdim2,fftable,nfftable,ffwork,nffwork) call pubz3d(isign,nfft1,nfft2,nfft3,array,#else $ ffwork,isys) $ nfftdim1,nfftdim2,array,nfftdim1,nfftdim2,fftable, call CCFFT3D(isign,nfft1,nfft2,nfft3,scale,array, isys(4)=0 isys(3)=0 isys(2)=0 isys(1)=3 scale = 1.d0#ifdef CRAY#endif $ nfftdim1,nfftdim2,fftable) call ZFFT3D(isign,nfft1,nfft2,nfft3,array,#ifdef SGIFFT isign = 1 integer nfftable,nffwork,isys(4) double precision scale integer isign,inc1,inc2,inc3 integer nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3 double precision array(*),fftable(*),ffwork(*) implicit none $ nfftable,nffwork) $ nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3, subroutine fft_forward(array,fftable,ffwork,c----------------------------------------------------------- end return#endif call pubz3di(nfft1,nfft2,nfft3,fftable,nfftable) write(6,*)'using public domain fft code'#else $ ffwork,isys) $ nfftdim1,nfftdim2,array,nfftdim1,nfftdim2,fftable, call CCFFT3D(isign,nfft1,nfft2,nfft3,scale,array, isys(4)=0 isys(3)=0 isys(2)=0 isys(1)=3 scale = 1.d0 isign = 0 write(6,*)'using cray fft code'#ifdef CRAY#endif call ZFFT3DI(nfft1,nfft2,nfft3,fftable)#ifdef SGIFFT double precision scale integer isign,inc1,inc2,inc3 integer nfftable,nffwork,isys(4) integer nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3 double precision array(*),fftable(*),ffwork(*) implicit none $ nfftable,nffwork) $ nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3, subroutine fft_setup(array,fftable,ffwork, end return#endif sizffwrk = 2*nfftmax sizfftab = 3*nfftable nffwork = nfftmax nfftable = 4*nfftmax + 15#else sizffwrk = nffwork sizfftab = nfftable nffwork = 4*nfftdim1*nfftdim2*nfftdim3 nfftable = 2*(nfftdim1+nfftdim2+nfftdim3+50)#ifdef CRAY#endif sizffwrk = nffwork sizfftab = nfftable nffwork = 0 nfftable = 2*(nfftdim1+nfftdim2+nfftdim3+50)#ifdef SGIFFT if ( nfft3 .eq. 2*n )nfftdim3 = nfft3+1 n = nfft3/2 nfftdim3 = nfft3 if ( nfft2 .eq. 2*n )nfftdim2 = nfft2+1 n = nfft2/2 nfftdim2 = nfft2 if ( nfft1 .eq. 2*n )nfftdim1 = nfft1+1 n = nfft1/2 nfftdim1 = nfft1 nfftmax = max(nfft1,nfft2,nfft3) integer n,nfftmax $ nfftable,nffwork,sizfftab,sizffwrk integer nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3, implicit none $ sizfftab,sizffwrk) $ nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork, subroutine get_fftdims(nfft1,nfft2,nfft3,C FFT CALLSc------------------------------------- end return10 continue d(j) = c(j-1) - c(j) do 10 j = 2,order d(1) = -c(1) integer j integer order double precision c(*),d(*) implicit none subroutine diff(c,d,order)c------------------------------------- end return c(1) = div*(1-x)*c(1)100 continue c(k-j) = div*((x+j)*c(k-j-1) + (k-j-x)*c(k-j)) do 100 j = 1,k-2 c(k) = div*x*c(k-1) div = 1.d0 / (k-1) integer j double precision div integer k double precision c(*),x implicit none subroutine one_pass(c,x,k)c------------------------------------- end return c(1) = 1.d0 - x c(2) = x c(order) = 0.d0 double precision c(order),x integer order implicit none subroutine init(c,x,order)c--------------------------------------------------- end return call one_pass(array,w,order)c one more recursion call diff(array,darray,order)c perform standard b-spline differentiation10 continue call one_pass(array,w,k) do 10 k = 3,order-1c compute standard b-spline recursion call init(array,w,order)c do linear case integer k double precision w,array(order),darray(order) integer order implicit nonec---------- use standard B-spline recursions: see doc file subroutine fill_bspline(w,order,array,darray)c--------------------------------------------------- end return100 continue call fill_bspline(w,order,theta3(1,n),dtheta3(1,n)) w = fr3(n)-int(fr3(n)) call fill_bspline(w,order,theta2(1,n),dtheta2(1,n)) w = fr2(n)-int(fr2(n)) call fill_bspline(w,order,theta1(1,n),dtheta1(1,n)) w = fr1(n)-int(fr1(n)) do 100 n = 1,numatoms integer n double precision w $ dtheta2(order,numatoms),dtheta3(order,numatoms) $ theta3(order,numatoms),dtheta1(order,numatoms), double precision theta1(order,numatoms),theta2(order,numatoms), double precision fr1(numatoms),fr2(numatoms),fr3(numatoms) integer numatoms,order implicit nonec---------------------------------------------------------------------c dtheta1,dtheta2,dtheta3: the 1st deriv of spline coeff arraysc theta1,theta2,theta3: the spline coeff arraysc OUTPUTc order: the order of spline interpolationc fr1,fr2,fr3 the scaled and shifted fractional coordsc numatoms: number of atomsc INPUT:c--------------------------------------------------------------------- $ theta1,theta2,theta3,dtheta1,dtheta2,dtheta3) $ numatoms,fr1,fr2,fr3,order, subroutine get_bspline_coeffs(c---------------------------------------------------------------------C SPLINEc---------------------------------------------------- end return300 continue290 continue280 continue w(i,j,k) = work(k) do 280 k = 1,n3 if ( isign .eq. 1) call cfftb(n3,work,table(1,3)) if ( isign .eq. -1) call cfftf(n3,work,table(1,3))270 continue work(k) = w(i,j,k) do 270 k = 1,n3 do 290 j = 1, n2 do 300 i = 1, n1cc transform along Z finally ...c200 continue190 continue180 continue w(i,j,k) = work(j) do 180 j = 1,n2 if ( isign .eq. 1) call cfftb(n2,work,table(1,2)) if ( isign .eq. -1) call cfftf(n2,work,table(1,2))170 continue work(j) = w(i,j,k) do 170 j = 1,n2 do 190 i = 1,n1 do 200 k = 1,n3cc transform along Y then ...c100 continue90 continue80 continue w(i,j,k) = work(i) do 80 i = 1,n1 if ( isign .eq. 1) call cfftb(n1,work,table(1,1)) if ( isign .eq. -1) call cfftf(n1,work,table(1,1))70 continue work(i) = w(i,j,k) do 70 i = 1,n1 do 90 j = 1, n2 do 100 k = 1, n3cc transform along X first ...cc nwork should be max(n1,n2,n3)c ntable should be 4*max(n1,n2,n3) +15 integer i,j,k double precision table(ntable,3) double complex work( nwork) double complex w(ld1,ld2,n3) integer n1,n2,n3,ld1,ld2,isign,ntable,nwork implicit none $ work,nwork) subroutine pubz3d(isign,n1,n2,n3,w,ld1,ld2,table,ntable,***************************************************************************** end return call cffti(n3,table(1,3)) call cffti(n2,table(1,2)) call cffti(n1,table(1,1))c ntable should be 4*max(n1,n2,n3) +15 double precision table(ntable,3) integer n1,n2,n3,ntable implicit none subroutine pubz3di(n1,n2,n3,table,ntable)c------------------------------------------------------------- end return write(6,*)'rms force err = ',rms rms = dsqrt(rms_num/rms_den)100 continue rms_den = rms_den + fdx(i)**2 + fdy(i)**2 + fdz(i)**2 $ (fz1(i)-fz2(i))**2 rms_num = rms_num + (fx1(i)-fx2(i))**2 + (fy1(i)-fy2(i))**2 + do 100 i = 1,numatoms rms_den = 0.d0 rms_num = 0.d0 integer i double precision rms_num,rms_den,rms $ fdx(*),fdy(*),fdz(*) double precision fx1(*),fy1(*),fz1(*),fx2(*),fy2(*),fz2(*), integer numatoms $ fx1,fy1,fz1,fx2,fy2,fz2,fdx,fdy,fdz) subroutine check_force(numatoms,c------------------------------------------------------------- end return write(6,*)'rel error = ',relerr write(6,*)'trace vir = ',svir write(6,*)'tot ene = ',etot relerr = 2.d0*abs(etot+svir)/(abs(etot)+abs(svir)) $ adj_vir(6)+rec_vir(6)+dir_vir(6) $ adj_vir(4)+rec_vir(4)+dir_vir(4)+ svir = adj_vir(1)+rec_vir(1)+dir_vir(1)+ etot = self_ene+adj_ene+dir_ene+rec_ene double precision etot,svir,relerr double precision adj_vir(6),rec_vir(6),dir_vir(6) double precision self_ene,adj_ene,dir_ene,rec_ene implicit none $ adj_vir,rec_vir,dir_vir) subroutine check_virial(self_ene,adj_ene,dir_ene,rec_ene,c------------------------------------------------------------- end return400 continue fz(n) = fz(n) + recip(3,1)*f1+recip(3,2)*f2+recip(3,3)*f3 fy(n) = fy(n) + recip(2,1)*f1+recip(2,2)*f2+recip(2,3)*f3 fx(n) = fx(n) + recip(1,1)*f1+recip(1,2)*f2+recip(1,3)*f3200 continue150 continue100 continue $ theta2(ith2,n) * dtheta3(ith3,n) f3 = f3 - nfft3 * term * theta1(ith1,n) * $ dtheta2(ith2,n) * theta3(ith3,n) f2 = f2 - nfft2 * term * theta1(ith1,n) * $ theta2(ith2,n) * theta3(ith3,n) f1 = f1 - nfft1 * term * dtheta1(ith1,n) *c force is negative of grad term = charge(n)*Q(1,i,j,k) i = i0 + 1 + (nfft1 - isign(nfft1,i0))/2 i0 = i0 + 1 do 100 ith1 = 1,order i0 = int(fr1(n)) - order j = j0 + 1 + (nfft2 - isign(nfft2,j0))/2 j0 = j0 + 1 do 150 ith2 = 1,order j0 = int(fr2(n)) - order k = k0 + 1 + (nfft3 - isign(nfft3,k0))/2 k0 = k0 + 1 do 200 ith3 = 1,order k0 = int(fr3(n)) - order f3 = 0.d0 f2 = 0.d0 f1 = 0.d0 do 400 n = 1,numatomsC$& nfft1,nfft2,nfft3,theta1,theta2,theta3,dtheta1,dtheta2,dtheta3)C$& SHARE(numatoms,fr1,fr2,fr3,charge,Q,fx,fy,fz,recip,order,C$DOACROSS LOCAL(f1,f2,f3,k0,k,j0,j,i0,i,term,n,ith1,ith2,ith3), double precision f1,f2,f3,term integer n,ith1,ith2,ith3,i0,j0,k0,i,j,k double precision Q(2,nfftdim1,nfftdim2,nfftdim3) $ dtheta3(order,numatoms) double precision dtheta1(order,numatoms),dtheta2(order,numatoms), $ theta3(order,numatoms),charge(numatoms) double precision theta1(order,numatoms),theta2(order,numatoms), double precision fx(numatoms),fy(numatoms),fz(numatoms) double precision fr1(numatoms),fr2(numatoms),fr3(numatoms) double precision recip(3,3) integer nfftdim1,nfftdim2,nfftdim3 integer numatoms,order,nfft1,nfft2,nfft3 implicit none $ order,nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3,Q) $ dtheta1,dtheta2,dtheta3,fx,fy,fz,fr1,fr2,fr3, $ numatoms,charge,recip,theta1,theta2,theta3, subroutine grad_sum(c----------------------------------------------------------- end return10 continue Q(i) = 0.d0 do 10 i = 1,ntot integer i double precision Q(ntot) integer ntot subroutine clearQ(Q,ntot)c----------------------------------------------------------- end return300 continue200 continue150 continue100 continue Q(1,i,j,k) = Q(1,i,j,k) + theta1(ith1,n) * prod i = i0 + 1 + (nfft1 - isign(nfft1,i0))/2 i0 = i0 + 1 do 100 ith1 = 1,order i0 = int(fr1(n)) - order prod = theta2(ith2,n)*theta3(ith3,n)*charge(n) j = j0 + 1 + (nfft2 - isign(nfft2,j0))/2 j0 = j0 + 1 do 150 ith2 = 1,order j0 = int(fr2(n)) - order k = k0 + 1 + (nfft3 - isign(nfft3,k0))/2 k0 = k0 + 1 do 200 ith3 = 1,order k0 = int(fr3(n)) - order do 300 n = 1,numatoms call clearQ(Q,ntot) ntot = 2*nfftdim1*nfftdim2*nfftdim3 double precision prod integer n,ntot,ith1,ith2,ith3,i0,j0,k0,i,j,k double precision Q(2,nfftdim1,nfftdim2,nfftdim3) $ theta3(order,numatoms),charge(numatoms) double precision theta1(order,numatoms),theta2(order,numatoms), double precision fr1(numatoms),fr2(numatoms),fr3(numatoms) integer nfftdim1,nfftdim2,nfftdim3 integer numatoms,order,nfft1,nfft2,nfft3 implicit nonec---------------------------------------------------------------------c Q the charge gridc OUTPUT:c order: the order of spline interpolationc nfftdim1,nfftdim2,nfftdim3: physical charge grid dimsc nfft1,nfft2,nfft3: the charge grid dimensionsc fr1,fr2,fr3 the scaled and shifted fractional coordsc theta1,theta2,theta3: the spline coeff arraysc charge: the array of atomic chargesc numatoms: number of atomsc INPUT:c--------------------------------------------------------------------- $ order,nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3,Q) $ numatoms,charge,theta1,theta2,theta3,fr1,fr2,fr3, subroutine fill_charge_grid(c------------------------------------------------------------------------ end return400 continue $ bsp_mod(k) = 0.5d0*(bsp_mod(k-1) + bsp_mod(k+1)) if ( bsp_mod(k) .lt. tiny ) do 400 k = 1,nfft300 continue bsp_mod(k) = sum1**2 + sum2**2250 continue sum2 = sum2 + bsp_arr(j)*dsin(arg) sum1 = sum1 + bsp_arr(j)*dcos(arg) arg = twopi*(k-1)*(j-1)/nfft do 250 j = 1,nfft sum2 = 0.d0 sum1 = 0.d0 do 300 k = 1,nfft tiny = 1.d-7 twopi = 2.d0*3.14159265358979323846 double precision sum1,sum2,twopi,arg,tiny integer j,kc storing it into bsp_modc Computes the modulus of the discrete fourier transform of bsp_arr, double precision bsp_mod(nfft),bsp_arr(nfft) integer nfft implicit none subroutine DFTMOD(bsp_mod,bsp_arr,nfft)c------------------------------------------------------------------------ end return call DFTMOD(bsp_mod3,bsp_arr,nfft3) call DFTMOD(bsp_mod2,bsp_arr,nfft2) call DFTMOD(bsp_mod1,bsp_arr,nfft1)150 continue bsp_arr(i) = array(i-1) do 150 i = 2,order+1100 continue bsp_arr(i) = 0.d0 do 100 i = 1,maxn call fill_bspline(w,order,array,darray) w = 0.d0 endif stop write(6,*)'nfft1-3 too large! check on MAXNFFT' if ( maxn .gt. MAXNFFT )then maxn = max(nfft1,nfft2,nfft3) endif stop write(6,*)'order too large! check on MAXORDER' if ( order .gt. MAXORDER )thenc Order is the order of the B spline approx.c bsp_mod1-3 hold these values, nfft1-3 are the grid dimensions,c this routine loads the moduli of the inverse DFT of the B splines integer i,maxn double precision bsp_arr(MAXNFFT) double precision array(MAXORDER),darray(MAXORDER),w parameter (MAXNFFT=1000) integer MAXNFFT parameter (MAXORDER=25) integer MAXORDER + bsp_mod3(nfft3) double precision bsp_mod1(nfft1),bsp_mod2(nfft2), integer nfft1,nfft2,nfft3,order implicit none $ nfft1,nfft2,nfft3,order) subroutine load_bsp_moduli(bsp_mod1,bsp_mod2,bsp_mod3,c--------------------------------------------------------------- end return100 continue fr3(n) = nfft3*(w - anint(w) + 0.5d0) w = x(n)*recip(1,3)+y(n)*recip(2,3)+z(n)*recip(3,3) fr2(n) = nfft2*(w - anint(w) + 0.5d0) w = x(n)*recip(1,2)+y(n)*recip(2,2)+z(n)*recip(3,2) fr1(n) = nfft1*(w - anint(w) + 0.5d0) w = x(n)*recip(1,1)+y(n)*recip(2,1)+z(n)*recip(3,1) do 100 n = 1,numatoms double precision w integer n double precision fr1(numatoms),fr2(numatoms),fr3(numatoms) double precision x(numatoms),y(numatoms),z(numatoms),recip(3,3) integer numatoms,nfft1,nfft2,nfft3c fr1,fr2,fr3 the scaled and shifted fractional coordsc OUTPUT:c recip: the 3x3 array of reciprocal vectors stored as columnsc x,y,z: arrays of cartesian coordsc numatoms: number of atomsc INPUT: implicit none $ fr1,fr2,fr3) $ numatoms,x,y,z,recip,nfft1,nfft2,nfft3, subroutine get_scaled_fractionals(c---------------------------------------------------------------------- end return $ order,nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3,Q) $ dtheta1,dtheta2,dtheta3,dx,dy,dz,fr1,fr2,fr3, $ numatoms,charge,recip,theta1,theta2,theta3, call grad_sum( $ nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork) $ Q,fftable,ffwork,nfft1,nfft2,nfft3, call fft_forward( & ) $ nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3,eer,virial,rkcut $ Q,ewald_coeff,volume,recip,bsp_mod1,bsp_mod2,bsp_mod3, call scalar_sum( $ nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork) $ Q,fftable,ffwork,nfft1,nfft2,nfft3, call fft_back( $ nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3,Q) $ numatoms,charge,theta1,theta2,theta3,fr1,fr2,fr3,order, call fill_charge_grid( $ theta1,theta2,theta3,dtheta1,dtheta2,dtheta3) $ numatoms,fr1,fr2,fr3,order, call get_bspline_coeffs( $ numatoms,x,y,z,recip,nfft1,nfft2,nfft3,fr1,fr2,fr3) call get_scaled_fractionals( $ nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork,sfft,sffw) call get_fftdims(nfft1,nfft2,nfft3,c get some integer array dimensions integer nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork,sfft,sffw $ fr2(numatoms),fr3(numatoms) $ dtheta2(siztheta),dtheta3(siztheta),fr1(numatoms), $ theta2(siztheta),theta3(siztheta),dtheta1(siztheta), double precision Q(siz_Q),ffwork(sizffwrk),theta1(siztheta),c STACK STORAGE: These arrays can be tossed after leaving this routine $ bsp_mod3(nfft3),fftable(sizfftab) double precision bsp_mod1(nfft1),bsp_mod2(nfft2),c HEAP STORAGE: These arrays need to be preserved throughout simulation integer sizfftab,sizffwrk,siztheta,siz_Qc SIZES of some arrays $ virial(3,3) double precision eer,dx(numatoms),dy(numatoms),dz(numatoms),c computed herec rigid molecule virial needs a correction term notc virial: virial due to k-space sum (valid for atomic scaling;c dx,dy,dz: forces incremented by k-space sumc eer: ewald reciprocal or k-space energyc OUTPUT $ charge(numatoms),recip(3,3),volume,ewald_coeff,rkcut double precision x(numatoms),y(numatoms),z(numatoms), integer numatoms,order,nfft1,nfft2,nfft3c nfft1,nfft2,nfft3: the dimensions of the charge grid arrayc and at least 4.c fifth degree is order 6 etc. The order must be an even number c order: the order of Bspline interpolation. E.g. cubic is order 4c ewald_coeff: ewald convergence parameterc volume: the volume of the unit cellc recip: 3x3 array of reciprocal unit cell vectors (stored as columns)c charge atomic chargesc x,y,z: atomic coordsc numatoms: number of atomsc INPUT implicit none $ fr1,fr2,fr3,rkcut) $ Q,ffwork,theta1,theta2,theta3,dtheta1,dtheta2,dtheta3, $ bsp_mod1,bsp_mod2,bsp_mod3,fftable, $ sizfftab,sizffwrk,siztheta,siz_Q, $ eer,dx,dy,dz,virial, $ order,nfft1,nfft2,nfft3, $ numatoms,x,y,z,charge,recip,volume,ewald_coeff, subroutine do_pmesh_kspace(c---------------------------------------------------- end return $ nfftable,nffwork) $ nfft1,nfft2,nfft3,nfftdim1,nfftdim2,nfftdim3, call fft_setup(dummy,fftable,ffwork, $ nfft1,nfft2,nfft3,order) call load_bsp_moduli(bsp_mod1,bsp_mod2,bsp_mod3, $ nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork,sfft,sffw) call get_fftdims(nfft1,nfft2,nfft3, integer nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork,sfft,sffw double precision dummy double precision fftable(sizfftab),ffwork(sizffwrk) + bsp_mod3(nfft3) double precision bsp_mod1(nfft1),bsp_mod2(nfft2), integer nfft1,nfft2,nfft3,order,sizfftab,sizffwrkc see DO_PMESH_KSPACE for explanation of arguments implicit none $ nfft1,nfft2,nfft3,order,sizfftab,sizffwrk) $ bsp_mod1,bsp_mod2,bsp_mod3,fftable,ffwork, subroutine pmesh_kspace_setup(c---------------------------------------------------- end returnC write(6,*)'total STACK storage needed = ',sizstackC write(6,*)'total HEAP storage needed = ',sizheap sizstack = siz_Q+6*siztheta+sizffwrk+3*numatoms sizheap = nfft1+nfft2+nfft3+sizfftab siz_Q = 2*nfftdim1*nfftdim2*nfftdim3 siztheta = numatoms*order $ sizfftab,sizffwrk) $ nfftdim1,nfftdim2,nfftdim3,nfftable,nffwork, call get_fftdims(nfft1,nfft2,nfft3, integer nfftdim1,nfftdim2,nfftdim3,nfftable,nffworkc heap or stack allocation.c This routine computes the above output parameters needed for c sizstack is total size of temporary storagec sizheap is total size of permanent storagec siztheta is size of arrays theta1-3 dtheta1-3c sizffwrk is temporary 3d fft work storagec sizfftab is permanent 3d fft table storagec sizfftab,sizffwrk,siztheta,siz_Qc OUTPUTc order is the order of B-spline interpolationc numatoms is number of atomsc nfft1,nfft2,nfft3 are the dimensions of the charge grid arrayc nfft1,nfft2,nfft3,numatoms,orderc INPUT $ sizfftab,sizffwrk,siztheta,siz_Q,sizheap,sizstack integer nfft1,nfft2,nfft3,numatoms,order, implicit none $ sizfftab,sizffwrk,siztheta,siz_Q,sizheap,sizstack) $ nfft1,nfft2,nfft3,numatoms,order, subroutine pmesh_kspace_get_sizes(c*-* c*-* Email:marchi@villon.saclay.cea.frc*-* Contact for info M. Marchi, CEA, Gif Sur Yvette 91191 (FRANCE) c*-* Full copyright notice at http://www.chim.unifi.it/orac/copyright4.0.htmlc*-* Copyright (C) 2000 Massimo Marchi and Piero Procaccic*-*
POR FAVOR LEA !! Hola chicos !!! Soy Caro, vivo en Ohio, EE. UU. Tengo 32 años, estoy muy feliz de haber recibido mi tarjeta de cajero automático en blanco de Adriano. Mi tarjeta de cajero automático en blanco puede retirar $ 4,000 por día. Lo obtuve de Él la semana pasada y ahora he retirado alrededor de $ 10,000 gratis. El cajero automático en blanco retira dinero de cualquier cajero automático y no tiene nombre porque está en blanco, solo su PIN estará en él, no se puede rastrear y ahora tengo dinero para negocios, compras y suficiente dinero para mí y mi familia. vivo. Estoy muy contento y feliz de haber conocido a Adriano porque conocí a cinco personas antes que él y no pudieron ayudarme. Pero estoy feliz ahora que Adriano envió la tarjeta a través de DHL y la recibí en dos días. Obtenga su propia tarjeta de él en este momento, la está dando por una pequeña tarifa para ayudar a las personas, incluso si es ilegal, pero ayuda mucho y nadie es atrapado o rastreado. Estoy feliz y agradecido con Adriano porque cambió mi historia de repente. La tarjeta funciona en todos los países. Es una buena noticia. La dirección de correo electrónico de Adriano es adrianohackers01@gmail.com
ResponderEliminar