335 lines
11 KiB
FortranFixed
335 lines
11 KiB
FortranFixed
program test
|
|
implicit none
|
|
integer ndim
|
|
double precision x0(10),xt(10),fvalue
|
|
|
|
ndim=2
|
|
x0(1)=-90.0d0
|
|
x0(2)=1000.0d0
|
|
call GSA(ndim,x0,xt,fvalue)
|
|
write(*,*)fvalue,xt(1),xt(2)
|
|
end
|
|
|
|
* Generalized Simulated Annealing - Code
|
|
* Program developed by Members of Lab. of Molecular Modeling
|
|
* Kleber C. Mundim (1995)
|
|
* ------------------------------------------------------
|
|
* Global optimization method using Generalized Simulated Annealing
|
|
* as for example;
|
|
* GSA Procedure + Your objective/coust function
|
|
* ______________ __________________
|
|
* | | Set of Parameters | |
|
|
* | | =================> | Routine with |
|
|
* | GSA-routine | | your objective |
|
|
* | | <================= | function |
|
|
* |______________| Objective function |__________________|
|
|
*
|
|
* ------------------------------------------------------
|
|
* First Version Jan./1994 (Kleber Mundim and Constantino Tsallis)
|
|
* Second Version Jun./1995 (Marcelo Moret)
|
|
* Third Version Sep./1995 (Thierry Lemaire and Amin Bassrei)
|
|
*
|
|
* Some Basic References and literature citation:
|
|
*
|
|
*1- Title : Geometry Optimization and Conformational Analysis Through
|
|
* Generalized Simulated Annealing
|
|
* Authors : Kleber C. Mundim and Constantino Tsallis
|
|
* Journal : Int.Journal of Quantum Chemistry, 58 (1996),373-381
|
|
*
|
|
*2- Title : Stochastic Molecular Optimization using Generalized
|
|
* Simulated Annealing
|
|
* Authors : Marcelo Moret, Pedro G. Pascutti, Paulo M. Bish
|
|
* and Kleber C. Mundim
|
|
* Journal : Journal of Computational Chemistry, 19 (1998) 647-657
|
|
*
|
|
*3- Title : Modeling Gravity Anomalies Through Generalized
|
|
* Simulated Annealing
|
|
* Authors : Kleber C. Mundim, Thierry Lemaire and Amin Bassrei
|
|
* Journal : Physica A, 252 (1998) 405-416
|
|
* ------------------------------------------------------
|
|
* Most important vectors and parameters used in the GSA routines.
|
|
* NDimension --> Number of parameters to be optimized (problem dependent).
|
|
* X_t(i) --> This vector contains the parameters.
|
|
* X_0(i) --> X_0(i) contains X_t at the time t-1.
|
|
* X_Min --> Minimal parameters obtained.
|
|
* To --> Initial Temperature.
|
|
* qV --> qV caracterize the Visiting Probability Function.
|
|
* qA --> qA Acceptance Parameter.
|
|
* qT --> qT Temperature parameter
|
|
* NStopMax --> Maximum number of GSA loops (problem dependent).
|
|
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
subroutine GSA(ndim,X0,Xt,fvalue)
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
* This routine start the GSA-loop
|
|
|
|
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
|
|
PARAMETER (MaxDim=500)
|
|
dimension X0(ndim),Xt(ndim)
|
|
|
|
COMMON /Par1/ qV1,qV2,qA1,qT1,D,exp1,exp2,Coef,Tqt,To,ToScale
|
|
COMMON /Par2/ qA,qV,qT
|
|
COMMON /Par3/ NDimension,NStopMax,NRAN
|
|
|
|
COMMON /Xvector/ X_t(Maxdim), X_0(Maxdim), X_Min(Maxdim)
|
|
|
|
DATA One /1.0D+00/
|
|
NDimension=ndim
|
|
CALL GSAini()
|
|
DO i=1,NDimension
|
|
X_0(i)=X0(i)
|
|
X_Min(i) = X_0(i)
|
|
X_t(i) = X_0(i)
|
|
ENDDO
|
|
func_0 = func(X_0)
|
|
func_Min = func_0
|
|
func_t = func_0
|
|
OneqA1 = One/qA1
|
|
Time = 0.0D0
|
|
NCycle = 0
|
|
DO WHILE (NCycle.LE.NStopMax)
|
|
Time = Time + One
|
|
NCycle = NCycle + 1
|
|
T = Tqt/((One+Time)**qT1 - One)
|
|
IF(D.EQ.0.0D0) THEN
|
|
Tup = One
|
|
ELSE
|
|
Tup = T**(D/(qT-3.0D0))
|
|
ENDIF
|
|
CALL Gfunc(T,Tup)
|
|
func_t = func(X_t)
|
|
IF(func_t .LE. func_0) THEN
|
|
DO I=1, NDimension
|
|
X_0(I) = X_t(I)
|
|
ENDDO
|
|
func_0 = func_t
|
|
IF(func_t .LE. func_Min) THEN
|
|
func_Min = func_t
|
|
DO I=1,NDimension
|
|
X_Min(I) = X_t(I)
|
|
ENDDO
|
|
ENDIF
|
|
ELSE
|
|
DeltaE = func_t - func_0
|
|
qA1= qA - One
|
|
PqA = One/((One+qA1*DeltaE/T)**OneqA1)
|
|
Rand = RAN3(NRAN)
|
|
IF(Rand .LT. PqA) THEN
|
|
DO I=1,NDimension
|
|
X_0(I) = X_t(I)
|
|
ENDDO
|
|
func_0 = func_t
|
|
ENDIF
|
|
ENDIF
|
|
ENDDO
|
|
fvalue=func_t
|
|
do I=1,ndim
|
|
Xt(I)=X_t(I)
|
|
enddo
|
|
return
|
|
END
|
|
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
SUBROUTINE Gfunc(T,Tup)
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
* This routine evaluate the new set o parameters using the
|
|
* visiting probability function g(q,T).
|
|
|
|
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
|
|
PARAMETER (MaxDim=500)
|
|
COMMON /Par1/ qV1,qV2,qA1,qT1,D,exp1,exp2,Coef,Tqt,To,ToScale
|
|
COMMON /Par3/ NDimension,NStopMax,NRAN
|
|
COMMON /Xvector/ X_t(Maxdim), X_0(Maxdim), X_Min(Maxdim)
|
|
|
|
DATA One /1.00000D+00/
|
|
|
|
DO I = 1,NDimension
|
|
R = RAN3(NRAN)
|
|
S = RAN3(NRAN)
|
|
DeltaX = Coef*Tup/(One+qV1*R*R/T**exp1)**exp2
|
|
IF(S.LE.0.5) DeltaX = -DeltaX
|
|
X_t(I) = X_0(I) + DeltaX
|
|
ENDDO
|
|
|
|
|
|
RETURN
|
|
END
|
|
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
FUNCTION dgamma(r)
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
* This routine evaluate usual Gamma function.
|
|
|
|
IMPLICIT REAL*8 (a - h, o - z)
|
|
|
|
PARAMETER (
|
|
1 p0 = 0.999999999999999990d+00, p1 = -0.422784335098466784d+00,
|
|
2 p2 = -0.233093736421782878d+00, p3 = 0.191091101387638410d+00,
|
|
3 p4 = -0.024552490005641278d+00, p5 = -0.017645244547851414d+00,
|
|
4 p6 = 0.008023273027855346d+00)
|
|
PARAMETER (
|
|
1 p7 = -0.000804329819255744d+00,p8 = -0.000360837876648255d+00,
|
|
2 p9 = 0.000145596568617526d+00,p10 = -0.000017545539395205d+00,
|
|
3 p11 = -0.000002591225267689d+00,p12 = 0.000001337767384067d+00,
|
|
4 p13 = -0.000000199542863674d+00)
|
|
|
|
n = NINT(r - 2)
|
|
w = r - (n + 2)
|
|
y = ((((((((((((p13 * w + p12)* w + p11)* w + p10)*
|
|
1 w + p9) * w + p8) * w + p7) * w + p6) * w + p5) *
|
|
2 w + p4) * w + p3) * w + p2) * w + p1) * w + p0
|
|
IF (n .gt. 0) THEN
|
|
w = r - 1
|
|
DO k = 2, n
|
|
w = w * (r - k)
|
|
END DO
|
|
ELSE
|
|
w = 1
|
|
DO k = 0, -n - 1
|
|
y = y * (r + k)
|
|
END DO
|
|
END IF
|
|
dgamma = w / y
|
|
|
|
RETURN
|
|
END
|
|
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
DOUBLE PRECISION FUNCTION RAN3(IDUM)
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
|
|
* This routine return an randomic number 0<= r <= 1
|
|
|
|
C IMPLICIT REAL*4(M)
|
|
C PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=2.5D-7)
|
|
PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1.D-9)
|
|
COMMON /RAN1A/ IFF,MJ,MK,INEXT,INEXTP,MA(55)
|
|
* DATA IFF /0/
|
|
* INEXT=0
|
|
* INEXTP=31
|
|
|
|
IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
|
|
IFF=1
|
|
MJ=MSEED-IABS(IDUM)
|
|
MJ=MOD(MJ,MBIG)
|
|
MA(55)=MJ
|
|
MK=1
|
|
DO 11 I=1,54
|
|
II=MOD(21*I,55)
|
|
MA(II)=MK
|
|
MK=MJ-MK
|
|
IF(MK.LT.MZ)MK=MK+MBIG
|
|
MJ=MA(II)
|
|
11 CONTINUE
|
|
DO 13 K=1,4
|
|
DO 12 I=1,55
|
|
MA(I)=MA(I)-MA(1+MOD(I+30,55))
|
|
IF(MA(I).LT.MZ)MA(I)=MA(I)+MBIG
|
|
12 CONTINUE
|
|
13 CONTINUE
|
|
INEXT=0
|
|
INEXTP=31
|
|
IDUM=1
|
|
ENDIF
|
|
|
|
INEXT=INEXT+1
|
|
IF(INEXT.EQ.56)INEXT=1
|
|
INEXTP=INEXTP+1
|
|
IF(INEXTP.EQ.56)INEXTP=1
|
|
MJ=MA(INEXT)-MA(INEXTP)
|
|
IF(MJ.LT.MZ)MJ=MJ+MBIG
|
|
MA(INEXT)=MJ
|
|
RAN3=MJ*FAC
|
|
|
|
RETURN
|
|
END
|
|
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
SUBROUTINE GSAini()
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
* This routine read and initialize the GSA parameters
|
|
* The number of parameters to be optimized (NDIMENSION) is problem
|
|
* dependent.
|
|
|
|
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
|
|
PARAMETER (MaxDim=500)
|
|
COMMON /Par1/ qV1,qV2,qA1,qT1,D,exp1,exp2,Coef,Tqt,To,ToScale
|
|
COMMON /Par2/ qA,qV,qT
|
|
COMMON /Par3/ NDimension,NStopMax,NRAN
|
|
|
|
DATA One /1.00000D+00/
|
|
CHARACTER tempo*8
|
|
|
|
qA=1.5d0
|
|
qT=1.5d0
|
|
qV=1.5d0
|
|
NStopMax=1000
|
|
To=1.0E-00
|
|
CALL TIME(tempo)
|
|
READ(tempo(7:8),"(I2)")NRAN
|
|
NRAN = -NRAN
|
|
|
|
D = 10.5d0
|
|
|
|
Pi = 3.14159265359D0
|
|
|
|
* Acceptance probability
|
|
qA1= qA - One
|
|
|
|
* Temperature
|
|
qT1 = qT - One
|
|
Tqt = To*(2.0D0**qT1-One)
|
|
|
|
* Visiting probability
|
|
qV1= qV - One
|
|
qV2= 2.0D0**qV1 - One
|
|
tmp= One/qV1 - 0.5D0
|
|
GamaDown = dgamma(tmp)
|
|
exp1 = 2.0D0/(3.0D0 - qV)
|
|
|
|
IF(D.EQ.0.0D0) THEN
|
|
Coef1 = One
|
|
exp2 = One/qv1 - 0.5D0
|
|
GamaUp = GamaDown
|
|
ELSE
|
|
Coef1 = (qV1/Pi)**(D*0.5D0)
|
|
exp2 = One/qV1 + 0.5D0*D - 0.5D0
|
|
GamaUp = dgamma(exp2)
|
|
ENDIF
|
|
|
|
Coef = Coef1*GamaUp/GamaDown
|
|
|
|
|
|
RETURN
|
|
END
|
|
|
|
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
FUNCTION func(X)
|
|
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * |
|
|
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
|
|
PARAMETER (MaxDim=500)
|
|
* This subrotuine make link between GSA and your Objective function
|
|
* X(nDimension) is a vector that contains the parameters set.
|
|
* GSA Procedure + Your objective/coust function
|
|
* ______________ __________________
|
|
* | | Set of Parameters X(i) | |
|
|
* | | =================> | Routine with |
|
|
* | GSA-routine | | your objective |
|
|
* | | <================= | function |
|
|
* |______________| Value of the function (f) |__________________|
|
|
|
|
COMMON /Par3/ NDimension,NStopMax,NRAN
|
|
|
|
DIMENSION X(MaxDim)
|
|
|
|
* CALL here your routine with the objective function
|
|
* X(i) is parameters set and f is the value of the objective function
|
|
! CALL funct(X,f)
|
|
! func = f
|
|
|
|
func=(X(1)-2.0d0)*(X(1)-2.0d0)+(X(2)-13.8d0)*(X(2)-13.8d0)
|
|
RETURN
|
|
END
|
|
|