我正在使用 PGI visual fortran 编写一个 fortran 程序。该代码由相当长的两个模块和一个主程序组成。对我来说似乎没有任何问题,程序编译和构建成功,没有警告或错误。但是在运行时或调试模式下会发生错误:
Signalled NONCONTINUABLE_EXCEPTION at 0x77d340f2 in function
RtlUnhandledExceptionFilter.
Cannot debug target program. This exception may be due to a
missing DLL required by the target program.
在此错误之后,程序的执行将停止。我尝试了所有我能猜到的方法来克服这个问题。但是,问题仍然存在。变量似乎定义正确。在调试模式下,所有变量似乎都正常。虽然变量 x、y 和 z 将由模块中的函数正确计算,但它们不能在主程序中计算,如下代码所示(通过文件 mainser.f90 中程序 main 的第 14 行的断点获得):
x: Children Could not be evalauted
我无法理解错误或访问冲突的来源。这是一个我不熟悉的已知问题吗?或者它是一个错误?我会感谢你的任何帮助。先感谢您。附上fortran语言代码。mainser.f90:
Program mainser
Use rand
Use ChainConfModule
Implicit NONE
Integer :: i, j, N, m, Nx, Ny, Nz, Nbox, box
Double Precision :: pi, eta, D, x(N,m), y(N,m), z(N,m)
! Definition of Constants
N=8;
m=6; ! Nx*Ny*Nz*4 must be equal to N*m
Nx=2; Ny=2; Nz=3;
Nbox=125; ! 1000 chains
pi=3.141592653589793D0
eta=pi/6.0D0*0.8D0
Call ChainConfig(N,m,Nx,Ny,Nz,eta,D,x,y,z)
End Program mainser
ChainConfModule.f90:
Module ChainConfModule
Use rand
Implicit None
Contains
!********************************************************************
Subroutine ChainConfig(N,m,Nx,Ny,Nz,eta,D,x,y,z)
Integer, Intent(in) :: N, m, Nx, Ny, Nz
Double Precision, Intent(in) :: eta
Double Precision, Intent(out):: x(N,m), y(N,m), z(N,m), D
Integer :: i,j,k,O
Double Precision :: max_eta, Xseg(N*m), Yseg(N*m), Zseg(N*m), Lx, Ly, Lz, D2, r2(N*m,N*m)
! Arranging Nm=N*m segments in a close-packed FCC stucture
max_eta=0.7404D0
i=FCCconf(N*m,Nx,Ny,Nz,max_eta,Xseg,Yseg,Zseg,Lx,Ly,Lz,D2,r2)
D=D2**(1.0D0/2.0D0);
If (i==0) Then
x=0; y=0; z=0;
return
End If
! Generating random config of PN chains (N chains each having m segments)
i=ChainGrowth(N,m,Xseg,Yseg,Zseg,r2,x,y,z)
If (i==0) Then
x=0; y=0; z=0;
return
End If
! Expand the system to desired density
Call expand(N,m,Lx,Ly,Lz,D,eta,x,y,z)
! Final Check for possible overlap of molecules
k=1;
Do i=1 , N
Do j=1 , m
Xseg(k)=x(i,j); Yseg(k)=y(i,j); Zseg(k)=z(i,j);
k=k+1;
End Do
End Do
O=checkoverlapall(Xseg,Yseg,Zseg,N*m,D2,r2);
If (O==1) Then
x=0; y=0; z=0;
return
End If
End Subroutine ChainConfig
!********************************************************************
Integer Function ChainGrowth(N,m,Xseg,Yseg,Zseg,r2,x,y,z)
Integer, intent(in) :: N, m
Double Precision, intent(in) :: Xseg(N*m), Yseg(N*m), Zseg(N*m)
Double Precision, intent(in) :: r2(N*m,N*m)
Double Precision, intent(out) :: x(N,m), y(N,m), z(N,m)
Integer :: i,j,Nm, neigh(N*m,N*m), XX(N*m,N*m), A(N*m), failed, iter, MaxIter, free(N*m), chain(N,m)
Integer :: temp(N), iter2, count, complete, temp2(m), temp3(N*m)
Double Precision :: minr
! Determine neighbor segments from their distances
Nm=N*m;
neigh=0;
minr=minval( r2(2:Nm,1) );
Do i=1 , Nm-1
Do j=i+1 , Nm
If (abs(r2(j,i)-minr)<1.0e-10) Then
neigh(i,j)=1;
neigh(j,i)=1;
End If
End Do
End Do
Do i=1,Nm
A(i)=i
End Do
Do i=1,Nm
XX(:,i)=A*neigh(:,i);
End Do
! Put Random Bonds between the segments
failed=1; iter=0;
maxiter=500000;
Do While (failed .AND. iter<maxiter)
! Select N segments randomly as heads of chains
free=1;
chain=0; temp=0;
Do i=1,Nm
j=0;
Do While (j==0 .OR. any(temp==j))
j=ceiling(grnd()*Nm);
End Do
temp(i)=j;
End Do
chain(:,1)=temp;
free(temp)=0;
! Select the sequence of next segments for all chains (independent growth of chains)
i=1; iter2=0; complete=0; count=0; A=0;
Do While (.NOT. complete .AND. iter2<100 )
j=chain(i,1)
temp3=free;
Call buildchain(j,N,m,XX,temp3,temp2)
chain(i,:)=temp2
If ( all(temp2 .NE. 0) .AND. different(A,temp2) ) Then
A(count+1:count+m)=chain(i,:)
count=count+m
free(chain(i,:))=0;
i=i+1;
If (i==N+1) Then
complete=1;
End If
Else
! Destroy the previous chain
If (i==1) Then
exit
Else
i=i-1;
free(chain(i,2:m))=1;
A(count-1:count-m)=0
count=count-m
iter2=iter2+1;
End If
End If
End Do
If (all(A .NE. 0)) Then
failed=0;
End If
iter=iter+1;
End Do ! while failed
! Construct x, y, z matrices
If (failed==0) Then
ChainGrowth=1;
Do i=1,N
Do j=1,m
x(i,j)=Xseg(chain(i,j));
y(i,j)=Yseg(chain(i,j));
z(i,j)=Zseg(chain(i,j));
End Do
End Do
Else
x=0; y=0; z=0;
ChainGrowth=0
End If
End Function ChainGrowth
!********************************************************************
Subroutine buildchain(i0,N,m,X,free,chain)
Integer, intent(in) :: i0, N, m, X(:,:)
Integer, intent(inout) :: free(:)
Integer, intent(out) :: chain(m)
Integer :: i, j, k, count, temp, c, alone
i=i0;
chain(1)=i0;
Do j=2,m
count=0;
alone=1;
Do c=1, N*m
temp=X(c,i)*free(c)
If (temp .NE. 0) Then
count=count+1;
If (free(temp)==1) Then
alone=0
End If
End If
End Do
If (alone==1) Then
chain(j:m)=0;
return
Else
k=ceiling(count*grnd())
count=0
Do c=1, N*m
temp=X(c,i)*free(c)
If (temp .NE. 0) Then
count=count+1;
If (count==k) Then
chain(j)=temp
i=temp
free(temp)=0
End If
End If
End Do
End If
End Do
End Subroutine buildchain
!********************************************************************
Integer function different(a,b)
Integer, intent(in) :: a(:),b(:)
Integer :: m,n,i,j
different=1;
m=size(a);
n=size(b);
Do i=1,m
Do j=1,n
If (a(i)==b(j)) Then
different=0;
return;
End If
End Do
End Do
End Function different
!********************************************************************
Integer Function FCCconf(N,Nx,Ny,Nz,eta,x,y,z,Lx,Ly,Lz,D2,r2)
Integer, Intent(in) :: N, Nx, Ny, Nz
Double Precision, Intent(in) :: eta
Double Precision, Intent(out):: x(N), y(N), z(N), Lx, Ly, Lz, D2, r2(N,N)
Integer :: i, cx, cy, cz, O
Double Precision :: pi, max_eta, a, V, D3, D
FCCconf=1
! Check the packing fraction to be under the maximum allowed
pi=3.141592653589793D0
max_eta=sqrt(2.0D0)*pi/6.0D0
If (eta>max_eta) Then
FCCconf=0
return;
End If
! Calculation of No. of unit cells in the box
If (N .NE. Nx*Ny*Nz*4) Then
FCCconf=0
return;
End If
! Filling the box with spheres in the FCC arrangement
a=1.0D0/DBLE(min(Nx, Ny, Nz)); ! length of a side of one cell
Lx=a*DBLE(Nx); Ly=a*DBLE(Ny); Lz=a*DBLE(Nz);
V=DBLE(Nx*Ny*Nz)*a**3.0D0;
D3=6.0D0*eta*V/(pi*DBLE(N));
D2=D3**(2.0D0/3.0D0);
D=D2**(1.0D0/2.0D0)
i=1;
Do cz=0 , Nz-1
Do cy=0 , Ny-1
Do cx=0 , Nx-1
x(i)=cx*a+0; y(i)=cy*a+0; z(i)=cz*a+0;
x(i+1)=cx*a+a/2; y(i+1)=cy*a+0; z(i+1)=cz*a+a/2;
x(i+2)=cx*a+0; y(i+2)=cy*a+a/2; z(i+2)=cz*a+a/2;
x(i+3)=cx*a+a/2; y(i+3)=cy*a+a/2; z(i+3)=cz*a+0;
i=i+4;
End Do
End Do
End Do
x=x+D/2.0D0; y=y+D/2.0D0; z=z+D/2.0D0
O=checkoverlapall(x,y,z,N,D2,r2);
If (O==1) Then
FCCconf=0
return;
End If
End Function FCCconf
!********************************************************************
Integer Function checkoverlapall(x,y,z,N,D2,r2)
! Check the overlap of all spheres in the box
Integer, intent(in) :: N
Double Precision, intent(in) :: D2, x(N), y(N), z(N)
Double Precision, intent(out) :: r2(N,N)
integer :: i, j
Double Precision :: dx, dy, dz
checkoverlapall=0;
! Check overlap of the generated configuration
Do i=1 , N-1
Do j=i+1 , N
dx=x(i)-x(j);
dy=y(i)-y(j);
dz=z(i)-z(j);
! NO minimum image convention
r2(j,i)=dx**2.0D0+dy**2.0D0+dz**2.0D0;
If (r2(j,i)-D2 < -1e-10) Then
checkoverlapall=1;
End If
End Do
End Do
End Function checkoverlapall
!********************************************************************
Subroutine expand(N,m,Lx,Ly,Lz,D,eta,x,y,z)
! Expand the close-packed structure to the desired density using a MC algorithm
Integer, intent(in) :: N,m
Double Precision, intent(in) :: D, eta
Double Precision, intent(inout) :: Lx,Ly,Lz, x(N,m),y(N,m),z(N,m)
Integer :: i, acc, accpt, cycl, Nc
Double Precision :: r, d_max, X_CM(N), Y_CM(N), Z_CM(N), X_CM_new, Y_CM_new, Z_CM_new
r=(0.74D0/eta)**(1.0D0/3.0D0);
Lx=r*Lx;
Ly=r*Ly;
Lz=r*Lz;
! Calculate center of mass of molecules
Do i=1,N
X_CM(i)=sum(x(i,:))/DBLE(m);
Y_CM(i)=sum(y(i,:))/DBLE(m);
Z_CM(i)=sum(z(i,:))/DBLE(m);
End Do
! MC algorithm for moving chains
Nc=1e6; ! No. of cycles
d_max=dble(D)/10.0D0;
Do cycl=1,Nc
accpt=0;
Do i=1,N
X_CM_new=X_CM(i)+(2*grnd()-1)*d_max;
Y_CM_new=Y_CM(i)+(2*grnd()-1)*d_max;
Z_CM_new=Z_CM(i)+(2*grnd()-1)*d_max;
acc=checkoverlapchain(N,m,i,X_CM_new,Y_CM_new,Z_CM_new,X_CM,Y_CM,Z_CM,x,y,z,D**2.0D0,Lx,Ly,Lz)
If (acc==1) Then
accpt=accpt+1;
End If
End Do
If (accpt/DBLE(N) < 0.5 ) Then
d_max=d_max*0.95;
Else
d_max=d_max*1.05;
End If
End Do
End Subroutine expand
!********************************************************************
Integer Function checkoverlapchain(N,m,index,X_CM_new,Y_CM_new,Z_CM_new,X_CM,Y_CM,Z_CM,x,y,z,D2,Lx,Ly,Lz)
! chain index moved with new CM coordinates
Integer, intent(in) :: N,m,index
Double Precision, intent(in) :: Lx, Ly, Lz, D2, X_CM_new, Y_CM_new, Z_CM_new
Double Precision, intent(inout) :: X_CM(N),Y_CM(N),Z_CM(N),x(N,m),y(N,m),z(N,m)
Integer :: i, j, k, O
Double Precision :: x_new(m),y_new(m),z_new(m), Xseg(N*m), Yseg(N*m), Zseg(N*m), r2(N*m,N*m), D
! Calculate new coordinates of chain segments
Do j=1,m
x_new(j)=x(index,j)-X_CM(index)+X_CM_new;
y_new(j)=y(index,j)-Y_CM(index)+Y_CM_new;
z_new(j)=z(index,j)-Z_CM(index)+Z_CM_new;
End Do
If ( any(x_new>Lx-D/2.0D0) .OR. any(x_new<D/2.0D0) .OR. any(y_new>Ly-D/2.0D0) .OR. any(y_new<D/2.0D0) &
.OR. any(z_new>Lz-D/2.0D0) .OR. any(z_new<D/2.0D0) ) Then
! Reject move
checkoverlapchain=0;
return
End If
k=1;
Do i=1 , N
If (i==index) Then
Do j=1 , m
Xseg(k)=x_new(j); Yseg(k)=y_new(j); Zseg(k)=z_new(j);
k=k+1;
End Do
Else
Do j=1 , m
Xseg(k)=x(i,j); Yseg(k)=y(i,j); Zseg(k)=z(i,j);
k=k+1;
End Do
End If
End Do
O=checkoverlapall(Xseg,Yseg,Zseg,N*m,D2,r2);
If (O==1) Then
! Reject move
checkoverlapchain=0;
Else
! Accept move
checkoverlapchain=1;
X_CM(index)=X_CM_new; Y_CM(index)=Y_CM_new; Z_CM(index)=Z_CM_new;
x(index,:)=x_new; y(index,:)=y_new; z(index,:)=z_new;
End If
End Function checkoverlapchain
!********************************************************************
End Module ChainConfModule
rand.f90:
! A C-program for MT19937: Real number version
! genrand() generates one pseudorandom real number (double)
! which is uniformly distributed on [0,1]-interval, for each
! call. sgenrand(seed) set initial values to the working area
! of 624 words. Before genrand(), sgenrand(seed) must be
! called once. (seed is any 32-bit integer except for 0).
! Integer generator is obtained by modifying two lines.
! Coded by Takuji Nishimura, considering the suggestions by
! Topher Cooper and Marc Rieffel in July-Aug. 1997.
!
! This library is free software; you can redistribute it and/or
! modify it under the terms of the GNU Library General Public
! License as published by the Free Software Foundation; either
! version 2 of the License, or (at your option) any later
! version.
! This library is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
! See the GNU Library General Public License for more details.
! You should have received a copy of the GNU Library General
! Public License along with this library; if not, write to the
! Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
! 02111-1307 USA
!
! Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
! When you use this, send an email to: matumoto@math.keio.ac.jp
! with an appropriate reference to your work.
!
!***********************************************************************
! Fortran translation by Hiroshi Takano. Jan. 13, 1999.
!
! genrand() -> double precision function grnd()
! sgenrand(seed) -> subroutine sgrnd(seed)
! integer seed
!
! This program uses the following non-standard intrinsics.
! ishft(i,n): If n>0, shifts bits in i by n positions to left.
! If n<0, shifts bits in i by n positions to right.
! iand (i,j): Performs logical AND on corresponding bits of i and j.
! ior (i,j): Performs inclusive OR on corresponding bits of i and j.
! ieor (i,j): Performs exclusive OR on corresponding bits of i and j.
!
!***********************************************************************
module rand
contains
subroutine sgrnd(seed)
implicit integer(a-z)
parameter(N = 624)
dimension mt(0:N-1)
! the array for the state vector
common /block/mti,mt
save /block/
! setting initial seeds to mt[N] using
! the generator Line 25 of Table 1 in
! [KNUTH 1981, The Art of Computer Programming
! Vol. 2 (2nd Ed.), pp102]
mt(0)= iand(seed,-1)
do 1000 mti=1,N-1
mt(mti) = iand(69069 * mt(mti-1),-1)
1000 continue
return
end subroutine sgrnd
!***********************************************************************
double precision function grnd()
implicit integer(a-z)
parameter(N = 624)
parameter(N1 = N+1)
parameter(M = 397)
parameter(MATA = -1727483681)
! constant vector a
parameter(UMASK = -2147483648)
! most significant w-r bits
parameter(LMASK = 2147483647)
! least significant r bits
! Tempering parameters
parameter(TMASKB= -1658038656)
parameter(TMASKC= -272236544)
dimension mt(0:N-1)
! the array for the state vector
common /block/mti,mt
save /block/
data mti/N1/
! mti==N+1 means mt[N] is not initialized
dimension mag01(0:1)
data mag01/0, MATA/
save mag01
! mag01(x) = x * MATA for x=0,1
TSHFTU(y)=ishft(y,-11)
TSHFTS(y)=ishft(y,7)
TSHFTT(y)=ishft(y,15)
TSHFTL(y)=ishft(y,-18)
if(mti.ge.N) then
! generate N words at one time
if(mti.eq.N+1) then
! if sgrnd() has not been called,
call sgrnd(4357)
! a default initial seed is used
endif
do 1000 kk=0,N-M-1
y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1)))
1000 continue
do 1100 kk=N-M,N-2
y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1)))
1100 continue
y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK))
mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1)))
mti = 0
endif
y=mt(mti)
mti=mti+1
y=ieor(y,TSHFTU(y))
y=ieor(y,iand(TSHFTS(y),TMASKB))
y=ieor(y,iand(TSHFTT(y),TMASKC))
y=ieor(y,TSHFTL(y))
if(y.lt.0) then
grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0)
else
grnd=dble(y)/(2.0d0**32-1.0d0)
endif
return
end function grnd
end module rand