store.f90

蒙特卡罗的一个程序分析 与大家分享 共同研究

subroutine StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
				      Conf, DXYZ, DROT, MaxMol, Ubins, Umin, Umax, &
					  tag_val, tag_mol, MaxEE, EESTEPS, EEW, &
				      LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
					  MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, &
					  Nmin, Nmax, Nham, LNW, Seed )

! This subroutine stores the current conformation to FinalConf if Istore = 0
! If Istore = 1, 2, or 3 the configuration is stored in a temporary file.

implicit none

integer, intent(inout)									:: Istore  

character*30, intent(in)							    :: FinalConf
	! If Istore==1,2 or 3 then the configuration is stored to a 
	! temporary file called 'tmp0'//'1,2, or 3'//'_'//TRIM(FinalConf).

integer, intent(in)										:: Nsp
integer, intent(in)										:: MaxSp

integer, dimension(0:MaxSp), intent(in)					:: Nmol

character*15, dimension(MaxSp), intent(in)				:: NAMEsp

integer, intent(in)										:: Conf

real, dimension(MaxSp), intent(in)						:: DXYZ
real, dimension(MaxSp), intent(in)						:: DROT

integer, intent(in)										:: MaxMol

integer, intent(in)										:: Ubins

real, intent(in)										:: Umin, Umax

integer, dimension(MaxSp), intent(in)					:: tag_val, tag_mol

integer, intent(in)										:: MaxEE
integer, dimension(MaxSp), intent(in)					:: EESTEPS
real, dimension(MaxEE,MaxSp), intent(in)				:: EEW

integer, dimension(MaxMol), intent(in)					:: LENGTHlj
integer, dimension(MaxMol), intent(in)					:: LENGTHion
integer, dimension(MaxMol), intent(in)					:: SPECIES

integer, intent(in)										:: MaxBeads

character*5, dimension(MaxBeads,MaxSp), intent(in)		:: BEADTYPE

integer, intent(in)										:: MaxLJ
integer, intent(in)										:: MaxIon

real, dimension(MaxLJ), intent(in)						:: Xlj, Ylj, Zlj
real, dimension(MaxIon), intent(in)						:: Xion, Yion, Zion

integer, intent(in)										:: Nmin, Nmax

integer, intent(in)										:: Nham

real, dimension(Nham), intent(in)						:: LNW

integer, intent(in)										:: Seed


! Local variables.

integer					:: i, j, k
integer					:: lenlj, lenion, sp
integer					:: ljb, ionb

integer, parameter		:: unit = 15

character*6, parameter  :: name01 = 'tmp01_'
character*6, parameter  :: name02 = 'tmp02_'
character*6, parameter  :: name03 = 'tmp03_'


if ( Istore == 0 ) then

     open( unit, file = FinalConf )

else if ( 1 <= Istore   .AND.   Istore <= 3 ) then

     ! Istore == 1 or 2 cover the usual case of saving the current
     ! configuration in two files, so that we have protection against
     ! a system crash.
     
	 if( Istore == 1 ) then
       
	    Istore = 2	    
        open( unit, file = name01//TRIM(FinalConf) )
     
	 else if ( Istore == 2 ) then
     
	    Istore = 3
        open( unit, file = name02//TRIM(FinalConf) )
     
	 else if ( Istore == 3 ) then
     
	    Istore = 1
        open( unit, file = name03//TRIM(FinalConf) )
     
	 end if

end if

! Write name and molecules of each species.
do i = 1, Nsp

     write(unit,"(A, T45, 4(I5,2X) )") NAMEsp(i), Nmol(i)

end do

! Write number of molecules per phase, BOXLENGTHs etc.
write(unit,*)
write(unit,*) Conf, Seed
write(unit,*)
write(unit,*) tag_val(1:Nsp), tag_mol(1:Nsp)
write(unit,*)

write(unit, "(A)") ' DXYZ, DROT '
    
write(unit,"( 10(1X,G24.16) )") DXYZ(1:Nsp)
    
write(unit,"( 10(1X,G24.16) )") DROT(1:Nsp)

write(unit,*)

write(unit, "(A)") '   Ubins        Umin              Umax '

write(unit,"( I8, 2(2x, G16.8))") Ubins, Umin, Umax

write(unit, "(A)") '    Nmin    Nmax '

write(unit,"(2I8)") Nmin, Nmax
	
write(unit,*)

write(unit,"(A)") ' HSGCMC Weights '

do i = 1, Nham

	write(unit,"(1x,G16.8E3)") exp( LNW(i) )
	
end do

!20 format( <Nham>(1x,G16.8E3) )

write(unit,*)

do i = 1, Nsp

	write(unit,"(A,I2)") ' EE Weights for Species ', i

	do j = 1, EESTEPS(i)

		write(unit,"(1x,G16.8E3)") EEW(j,i)
		
	end do

	write(unit,*)

end do

! Write configuration.

ljb = 0 
ionb = 0
     
 do j = 1, Nmol(0)

	lenlj = LENGTHlj(j)
	lenion = LENGTHion(j)

	sp = SPECIES(j)

	do k = 1, lenlj + lenion

		if( BEADTYPE(k,sp) == 'LJ' ) then

			ljb = ljb + 1

			write(unit, "(3(1X,G24.16), 3(1X,I5))")	&
				Xlj(ljb), Ylj(ljb), Zlj(ljb), j, k, sp
					
		else if( BEADTYPE(k,sp) == 'ION' ) then

			ionb = ionb + 1

			write(unit, "(3(1X,G24.16), 3(1X,I5))")	&
				Xion(ionb), Yion(ionb), Zion(ionb), j, k, sp

		end if

	end do

end do

close( unit )
	 
return

end subroutine StoreConf