8 integer*4 :: i,j,m,n,iloc,pos,iter,maxiter,numarg,tmplen
9 integer*4 :: WybierzOsobnika,Najlepszy,ZnajdzPodobnego,FindWorst
11 integer*4 :: mnoznik = 10
12 character*16 :: mode,argn,argi,dir_nam
13 character*15 :: filename
14 real*8 :: loc,fc,sumfitness,sump,avrp,maxfc
17 real*8 :: FunkcjaCeluB, FunkcjaCeluB2, FunkcjaCeluR
18 character(80),dimension(5) :: wagi
19 logical :: searching = .false.
20 logical :: logdata = .true.
21 logical :: fileex = .false.
22 logical :: retry = .false.
23 logical :: dir_ex = .false.
24 logical :: debug = .true.
25 character*200 :: text, tmptext,tmptext2
26 character(len=*), parameter :: FMT = "(19F10.5,E15.7,F10.5)"
28 c =======================================================================
30 c =======================================================================
33 call ReadInput(status)
34 if ((maxminstep.gt.0).and.(generation.eq.0)) then
40 write(*,*) "There ware input file errors - check log."
45 r = rand(timeArray(1)*10000+timeArray(2)*100+timeArray(3))
50 c-------------------------------------------------------
51 c bank & populacja arrays store this
53 c 20 - energy obtained from zscore
55 c-------------------------------------------------------
56 allocate (bank(banksize,21))
57 allocate (populacja(BANK_MULTIPLIER*banksize,21))
58 allocate (pairs(BANK_MULTIPLIER*banksize))
62 do i=1,BANK_MULTIPLIER*banksize
75 inquire(FILE=ostatefn,EXIST=fileex)
76 c No. Prepere next generation
78 open(ostate,file=ostatefn)
79 read(ostate,'(I4)') generation
80 read(ostate,'(L2)') do_optima
81 read(ostate,'(L2)') do_ga
82 read(ostate,'(F7.5)') avrd
85 call write2log("Doing GA in this step")
87 if (.not.do_optima) then
89 call write2log("ZSCORE weights minimalization disabled for now")
90 generation=generation+1
92 write(tmptext,'(I4)') generation
93 call write2log("This is genaration "//tmptext)
94 call ReadOptimaW(BANK_MULTIPLIER*banksize,populacja)
96 c Yes. Generate random zero-population
98 call GenPopulation(BANK_MULTIPLIER*banksize,populacja)
99 write(tmptext,'(I4)') generation
100 call write2log("This is genaration "//tmptext)
105 c End of "fisrst time here?" code
109 c Do we actual use a Genetic Algorithm?
112 if (do_ga.eq..true.) then
114 c Yes. This is only done when doing second pass with ZSCORE (without weight minimalizaton)
115 c or with weight minimalization disabled (maxmin=0) just to obtain final score)
117 c \\---//================================================================
119 c "// Genetic Algorithm code starts here
121 c =//---\\===============================================================
123 call ReadZEnergy(BANK_MULTIPLIER*banksize,populacja)
124 do i=1,BANK_MULTIPLIER*banksize
125 c populacja(i,20)=FunkcjaCeluR(populacja(i,20))
126 write(*,*) "Ind ",i,populacja(i,20)
129 call write2log("==[ Doing GA now ]=============================")
132 c Populate bank depending on algorithm
136 call GetNBest(populacja,bank,banksize)
138 call GetNBest(populacja,bank,banksize)
140 call GetNBest(populacja,bank,banksize)
142 c --- debug begin ---
144 write(*,*) "generation: ",generation
145 write(*,*) "maxminstep: ",maxminstep
146 write(*,*) "do_ga: ",do_ga
147 write(*,*) "do_optima: ",do_optima
148 write(*,*) "do_fs: ",do_fs
149 c write(*,*) "cicutoff: ", cicutoff
156 c Fill the bank just after the first time we get the score
158 if (((generation.eq.1).and.(maxminstep.eq.0)).or.((generation.eq&
159 &.1).and.(maxminstep.gt.0))) then
161 c --- debug begin ---
163 write(*,*) "Calling GetNBest..."
166 call GetNBest(populacja,bank,banksize)
167 call CalcAvgDist(bank,avrd)
168 write(tmptext,'(F7.5)') avrd
169 call write2log("Average distance between individuals in initial&
170 & bank is "//trim(tmptext))
171 csacutoff=(maxco*avrd)-generation*avrd*(maxco-minco)/maxgen
172 write(tmptext,'(F7.5)') csacutoff
173 call write2log("CSA cutoff is now set to "//trim(tmptext))
174 c csacutoff=maxco*avrd
178 c --- debug begin ---
181 write(*,FMT) bank(i,:)
186 call CalcFitness(bank,sumfitness)
194 c --- debug begin ---
196 write(*,*) "Szukam podobnego"
201 write(tmptext,'(F7.5)') avrd
202 call write2log("Average distance in bank is "//trim(tmptext))
204 do i=1,BANK_MULTIPLIER*banksize
205 write(tmptext,'(I4)') i
206 call write2log("Checking ind "//trim(tmptext))
207 j=ZnajdzPodobnego(banksize,bank,populacja(i,:),csacutoff)
209 if (populacja(i,20).lt.bank(j,20)) then
210 write(tmptext,'(I4)') j
211 write(tmptext2,'(I4)') i
212 call write2log("Swaping ind"//trim(tmptext)//" from bank to &
213 &ind "//trim(tmptext2)//" from population")
214 bank(j,:)=populacja(i,:)
217 j=FindWorst(banksize,bank)
218 write(tmptext,'(I4)') j
219 write(tmptext2,'(I4)') i
220 if (populacja(i,20).lt.bank(j,20)) then
221 call write2log("Worst in bank is "//trim(tmptext))
222 call write2log("Swaping worst ind in bank to "//trim(tmptext&
224 bank(j,:)=populacja(i,:)
229 c -----------------------------------------------------------------------
230 c Calculate fitness for reproduction
231 c -----------------------------------------------------------------------
232 call CalcFitness(bank,sumfitness)
234 c --- debug begin ---
236 write (*,*) "Dumping bank: "
238 write(*,FMT) bank(i,:)
240 write (*,*) "Sumfitness: ", sumfitness
246 csacutoff=maxco*avrd-generation*avrd*(maxco-minco)/maxgen
247 c csacutoff=csacutoff-(generation*cicutoff/maxgen)
248 c csacutoff=cicutoff*(0.8**(iter-1))
250 write(tmptext,'(F7.5)') csacutoff
251 call write2log("CSA cutoff is now set to "//trim(tmptext))
254 write(*,*) "Some stuff here in the future"
261 c Have we done last generation ?
262 if (generation.eq.maxgen) then
263 call write2log("Maximum number of genarations reached. QUITING.")
266 c Prapering for next generation so increase counter
269 c -----------------------------------------------------------------------
271 c -----------------------------------------------------------------------
272 call write2log("==[ Reproducting ]============================")
275 do i=1,BANK_MULTIPLIER*banksize
276 m=WybierzOsobnika(banksize,bank)
277 write(tmptext,'(I4)') m
278 write(tmptext2,'(I4)') i
279 call write2log("Reproducting individual "//trim(tmptext)//" from&
280 & bank as new individual "//trim(tmptext2))
281 populacja(i,:)=bank(m,:)
283 case('better','betterv2')
284 populacja(1,:)=bank(1,:)
285 call write2log("Reproducting first (best) individual")
286 do i=2,BANK_MULTIPLIER*banksize
287 m=WybierzOsobnika(banksize,bank)
288 write(tmptext,'(I4)') m
289 write(tmptext2,'(I4)') i
290 call write2log("Reproducting individual "//trim(tmptext)//" from&
291 & bank as new individual "//trim(tmptext2))
292 populacja(i,:)=bank(m,:)
296 c -----------------------------------------------------------------------
298 c -----------------------------------------------------------------------
299 do i=1,BANK_MULTIPLIER*banksize
303 do i=1,BANK_MULTIPLIER*banksize
304 if (pairs(i).eq.0) then
305 210 pos=1+int(rand(0)/(1/real(BANK_MULTIPLIER*banksize)))
310 if (pairs(j).eq.pos) then
318 call write2log("==[ Pairing ]==============================")
319 do i=1,BANK_MULTIPLIER*banksize
320 write(tmptext,'(A11,I3.3,A5,I3.3,A21)') "Individual ",i," and ",p&
321 &airs(i)," will be crossed over"
322 call write2log(trim(tmptext))
325 c -----------------------------------------------------------------------
327 c -----------------------------------------------------------------------
328 call write2log("==[ Crossing ]==============================")
329 do i=1,BANK_MULTIPLIER*banksize
330 if (pairs(i) .gt. i) then
331 call CrossoverS(populacja(i,:),populacja(pairs(i),:))
336 call write2log(" WLONG WSCP WELEC WBOND WANG &
337 &WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC WT&
338 &URN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR")
339 do i=1,BANK_MULTIPLIER*banksize
342 write(text,'(A4,I3.3,A2)') "Ind",i,": "
344 write(tmptext,'(F7.5)') populacja(i,j)
345 text = text(1:tmplen) // tmptext
348 call write2log(trim(text))
351 c -----------------------------------------------------------------------
353 c -----------------------------------------------------------------------
355 call write2log("==[ Mutating ]==============================")
357 if (IGNORE_WEIGHT(i).eq.1) then
358 write(tmptext,'(I2)') i
359 call write2log("Parameter "//trim(WNAME(i))//" will be ignored")
362 do i=1,BANK_MULTIPLIER*banksize
363 call MutationV(populacja(i,:),i,MUTRATIO)
365 c for betterv2 and csa restore best individual
367 case('betterv2','csa')
368 populacja(1,:)=bank(1,:)
372 c Genetic Algorithm blokcode stops here
376 c ---------------------------------
377 c Do some final stuff
378 c ---------------------------------
380 c Set the variables before writing to state file
382 if ((.not.do_ga).and.(.not.do_optima)) then
392 do_optima=.not.do_optima
393 if (generation.eq.0) then
406 c Write te current state and population
413 write(tmptext,'(I)') generation+1
414 call write2log("Preparing inputs for next generation ("//trim(adju&
416 call CreateInputs(BANK_MULTIPLIER*banksize,populacja)
418 c All done? Then let's quit.
424 c ======================================================================
426 c MAIN CODE ENDS HERE
428 c ======================================================================
432 2000 deallocate(populacja)
437 c ======================================================================
438 c Clustering subroutines
439 c ======================================================================
441 include 'cluster.inc'
444 c ======================================================================
445 c Write2log subroutine
446 c ======================================================================
448 c ----------------------------------------------------------------------
450 subroutine Write2log(text)
453 character*200 :: tmptext
454 character(len=*) :: text
456 call itime(timeArray)
457 inquire(file=ologfn,exist=ex)
459 open(olog, file=ologfn, access="append")
460 write(olog,"(A)") "==============================================&
462 tmptext = "= UNRES weights optimizer version "//version
463 write(olog,"(A)") tmptext
464 write(olog,"(A)") info
465 write(olog,"(A)") "==============================================&
469 open(olog, file = ologfn , access = "append")
470 write(olog,"(I2.2,A1,I2.2,A1,I2.2,A2,A)") timeArray(1),":",timeAr&
471 &ray(2),":",timeArray(3),": ",text
473 end subroutine Write2log
475 c ======================================================================
476 c GetNBest subroutine
477 c ======================================================================
478 c Fills up the bank population up to banksize with individuals from
479 c inputpop having the best(lowest) score
480 c ----------------------------------------------------------------------
482 subroutine GetNBest(inputpop,bank,banksize)
483 integer*4 :: banksize, i,j,idx
484 real*8 :: inputpop, bank, best, last
486 dimension inputpop(banksize*BANK_MULTIPLIER,21)
487 dimension bank(banksize,21)
493 do i=1,banksize*BANK_MULTIPLIER
494 if ((inputpop(i,20).lt.best).and.(inputpop(i,20).gt.last)) then
499 bank(j,:)=inputpop(idx,:)
502 end subroutine GetNBest
504 c ======================================================================
505 c FindWorst subroutine
506 c ======================================================================
507 c Finds the worst (highest score) individual in population
508 c ----------------------------------------------------------------------
510 integer*4 function FindWorst(rozmiar,pop)
511 integer*4 :: rozmiar, i, idx
513 dimension pop(rozmiar,21)
518 if (pop(i,20).gt.last) then
527 c ======================================================================
528 c GenPopulation subroutine
529 c ======================================================================
530 c Generates a population (pop) of nind individuals with 19 random
531 c wights from a set of ranges (see GA.inc).
532 c ----------------------------------------------------------------------
534 subroutine GenPopulation(nind,pop)
535 integer*4 :: i,j,nind
537 dimension pop(nind,21)
543 pop(j,1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
545 pop(j,2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
547 pop(j,3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
549 pop(j,4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
551 pop(j,5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
553 pop(j,6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
555 pop(j,7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
557 pop(j,8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
559 pop(j,9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
561 pop(j,10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
563 pop(j,11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
565 pop(j,12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
567 pop(j,13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
569 pop(j,14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
571 pop(j,15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
573 pop(j,16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
575 pop(j,17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
577 pop(j,18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
579 pop(j,19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
583 end subroutine GenPopulation
585 c ======================================================================
586 c CrossoverS subroutine
587 c ======================================================================
588 c Does a simple value crossing-over on a pair
589 c of two given individuals, giving two new individuals.
590 c ----------------------------------------------------------------------
592 subroutine CrossoverS(ind1,ind2)
593 real*8 :: ind1(19),ind2(19),temp(19)
595 loc = 1 + int(rand(0)/(1/18.0))
601 end subroutine CrossoverS
603 c ======================================================================
604 c MutationV subroutine
605 c ======================================================================
606 c Does a value mutation for the given individual (ind) with the
607 c given mutation ratio (mratio).
608 c ----------------------------------------------------------------------
610 subroutine MutationV(ind,id,mratio)
616 character*100 :: tmptext
617 character*150 :: logtext
620 if (IGNORE_WEIGHT(i).eq.1) then
621 c write(tmptext,'(I2)') i
622 c call write2log("Skipping weight "//trim(tmptext))
625 if (rand(0)<mratio) then
629 ind(1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
632 ind(2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
635 ind(3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
638 ind(4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
641 ind(5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
644 ind(6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
647 ind(7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
650 ind(8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
653 ind(9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
656 ind(10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
659 ind(11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
662 ind(12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
665 ind(13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
668 ind(14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
671 ind(15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
674 ind(16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
677 ind(17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
680 ind(18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
683 ind(19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
685 write(tmptext,'(I4,A13,I2,A7,F7.5,A5,F7.5)') id," at position "&
686 &,i," from: ",ti," to: ",ind(i)
687 logtext = "Mutation occured in individual "//tmptext
688 call write2log (logtext)
694 c ======================================================================
695 c FunkcjaCeluB function
696 c ======================================================================
697 c Simple (dummy) scoring function
698 c ----------------------------------------------------------------------
700 real*8 function FunkcjaCeluB(osobnik)
702 real*8 :: osobnik(20)
704 FunkcjaCeluB = 20 + -exp(-abs(WLONG-osobnik(1))) &
705 & - exp(-abs(WSCP-osobnik(2))) &
706 & - exp(-abs(WELEC-osobnik(3))) &
707 & - exp(-abs(WBOND-osobnik(4))) &
708 & - exp(-abs(WANG-osobnik(5))) &
709 & - exp(-abs(WSCLOC-osobnik(6))) &
710 & - exp(-abs(WTOR-osobnik(7))) &
711 & - exp(-abs(WTORD-osobnik(8))) &
712 & - exp(-abs(WCORRH-osobnik(9))) &
713 & - exp(-abs(WCORR4-osobnik(10))) &
714 & - exp(-abs(WCORR5-osobnik(11))) &
715 & - exp(-abs(WCORR6-osobnik(12))) &
716 & - exp(-abs(WEL_LOC-osobnik(13))) &
717 & - exp(-abs(WTURN3-osobnik(14))) &
718 & - exp(-abs(WTURN4-osobnik(15))) &
719 & - exp(-abs(WTURN6-osobnik(16))) &
720 & - exp(-abs(WVDWPP-osobnik(17))) &
721 & - exp(-abs(WHPB-osobnik(18))) &
722 & - exp(-abs(WSCCOR-osobnik(19)))
726 c ======================================================================
727 c FunkcjaCeluB2 function
728 c ======================================================================
729 c Simple (dummy) scoring function
730 c ----------------------------------------------------------------------
732 real*8 function FunkcjaCeluB2(osobnik)
734 real*8 :: osobnik(20)
736 FunkcjaCeluB2 = 1 + 4*(WLONG-osobnik(1))**2 &
737 & + 4*(WSCP-osobnik(2))**2 &
738 & + 4*(WELEC-osobnik(3))**2 &
739 & + 4*(WBOND-osobnik(4))**2 &
740 & + 4*(WANG-osobnik(5))**2 &
741 & + 4*(WSCLOC-osobnik(6))**2 &
742 & + 4*(WTOR-osobnik(7))**2 &
743 & + 4*(WTORD-osobnik(8))**2 &
744 & + 4*(WCORRH-osobnik(9))**2 &
745 & + 4*(WCORR4-osobnik(10))**2 &
746 & + 4*(WCORR5-osobnik(11))**2 &
747 & + 4*(WCORR6-osobnik(12))**2 &
748 & + 4*(WEL_LOC-osobnik(13))**2 &
749 & + 4*(WTURN3-osobnik(14))**2 &
750 & + 4*(WTURN4-osobnik(15))**2 &
751 & + 4*(WTURN6-osobnik(16))**2 &
752 & + 4*(WVDWPP-osobnik(17))**2 &
753 & + 4*(WHPB-osobnik(18))**2 &
754 & + 4*(WSCCOR-osobnik(19))**2
758 c ======================================================================
759 c Weights2str subroutine
760 c ======================================================================
761 c Writes weights of individual(osobnik) to string (lista)
762 c ----------------------------------------------------------------------
764 subroutine Weights2Str(osobnik,coff,lista)
765 real*8, dimension(19) :: osobnik
769 character(80) :: lista
777 write(a,"(A6,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)") "WLONG=",os&
778 &obnik(1)," WSCP=",osobnik(2)," WELEC=",osobnik(3)," WBOND=",osobni&
779 &k(4)," WANG=",osobnik(5)," &"
781 write(a,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)") "WSCLOC=",o&
782 &sobnik(6)," WTOR=",osobnik(7)," WTORD=",osobnik(8)," WCORRH=",osob&
783 &nik(9),"WCORR5=",osobnik(11)," &"
785 write(a,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)") "WCORR6=",o&
786 &sobnik(12)," WEL_LOC=",osobnik(13)," WTURN3=",osobnik(14)," WTURN4&
787 &=",osobnik(15),"WTURN6=",osobnik(16)," &"
789 write(a,"(A7,F7.5,A6,F7.5,A8,F7.5,A)") "WVDWPP=",osobnik(17)," WHP&
790 &B=",osobnik(18)," WSCCOR=",osobnik(19)," &
793 write(a,"(A7,F7.5,A8,F7.5)") "CUTOFF=",coff," WCORR4=",osobnik(10)
798 c ======================================================================
799 c FunkcjaCeluR function
800 c ======================================================================
801 c Real scoring function
802 c ----------------------------------------------------------------------
804 real*8 function FunkcjaCeluR(zscore)
806 if (zscore.eq.0) then
809 FunkcjaCeluR = 1/zscore
814 c ======================================================================
815 c ReadInput subroutine
816 c ======================================================================
818 subroutine ReadInput(status)
825 character*100 :: wiersz,tmp
826 inquire(FILE=inpfn,EXIST=ex)
830 call write2log('==[ Reading main input ]======================')
834 read(inp, '(A)', iostat=stat) wiersz
837 if ((wiersz(1:4).eq.'PDB=').or.(wiersz(1:4).eq.'pdb=')) then
839 tmp = wiersz(5:len_trim(wiersz))
841 if (tmp(i:i).eq.' ') then
845 if (npdb.gt.maxnpdb) then
846 call write2log("Number of input PDB exceeds maxnpdb!")
851 if (index(trim(tmp),' ').gt.0) then
852 pdbfiles(i)=tmp(1:index(trim(tmp),' '))
854 pdbfiles(i)=tmp(1:len_trim(tmp))
856 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
858 endif ! Koniec czytania "PDB="
861 if ((wiersz(1:4).eq.'ALG=').or.(wiersz(1:4).eq.'alg=')) then
862 alg=wiersz(5:len_trim(wiersz))
865 call write2log ("Using 'simple' algorithm")
867 call write2log ("Using 'better' algorithm")
869 call write2log ("Using 'betterv2' algorithm")
871 call write2log ("Using 'CSA' algorithm")
873 call write2log ("Using 'Clustering CSA' algorithm")
876 call write2log ("Unknown algorithm. Using 'csa' as default")
881 select case (wiersz(1:12))
882 case('GENERATIONS=','generations=')
883 tmp = wiersz(13:len_trim(wiersz))
884 read(tmp,'(I)') maxgen
888 c select case(wiersz(1:9))
889 c case('CICUTOFF=','cicutoff=')
890 c tmp = wiersz(10:len_trim(wiersz))
891 c read(tmp(1:len_trim(tmp)),'(F7.5)') cicutoff
892 c call write2log("Initial CSA cutoff is set to "//tmp)
896 select case(wiersz(1:6))
897 case('MINCO=','minco=')
898 tmp = wiersz(7:len_trim(wiersz))
899 read(tmp(1:len_trim(tmp)),'(F7.5)') minco
900 call write2log("Minimal CSA cutoff factor is set to "//tmp)
904 select case(wiersz(1:6))
905 case('MAXCO=','maxco=')
906 tmp = wiersz(7:len_trim(wiersz))
907 read(tmp(1:len_trim(tmp)),'(F7.5)') maxco
908 call write2log("Maximal CSA cutoff factor is set to "//tmp)
913 select case(wiersz(1:11))
914 case('POPULATION=','population=')
915 tmp = wiersz(12:len_trim(wiersz))
916 read(tmp,'(I)') banksize
917 call write2log("Bank size is set to "//tmp)
921 select case(wiersz(1:13))
922 case('WHAMTEMPLATE=','whamtemplate=')
923 tmp = wiersz(14:len_trim(wiersz))
926 if (tmp(i:i).eq.' ') then
930 if (ntwham.gt.maxnpdb) then
931 call write2log("Number of WHAM templates exceeds maxnpdb!")
936 if (index(trim(tmp),' ').gt.0) then
937 whamtemplate(i)=tmp(1:index(trim(tmp),' '))
939 whamtemplate(i)=tmp(1:len_trim(tmp))
941 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
946 select case(wiersz(1:14))
947 case('MREMDTEMPLATE=','mremdtemplate=')
948 tmp = wiersz(15:len_trim(wiersz))
951 if (tmp(i:i).eq.' ') then
955 if (ntmremd.gt.maxnpdb) then
956 call write2log("Number of MREMD templates exceeds maxnpdb!")
961 if (index(trim(tmp),' ').gt.0) then
962 mremdtemplate(i)=tmp(1:index(trim(tmp),' '))
964 mremdtemplate(i)=tmp(1:len_trim(tmp))
966 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
971 select case(wiersz(1:7))
972 case('MAXMIN=','maxmin=')
973 tmp = wiersz(8:len_trim(wiersz))
974 read(tmp,'(I)') maxminstep
975 if (maxminstep.gt.0) then
976 call write2log("ZSCORE weights minimalization set to "//tmp)
979 call write2log("ZSCORE weights minimalization disabled")
985 select case(wiersz(1:8))
986 case('SCRIPTS=','scripts=')
988 tmp = wiersz(9:len_trim(wiersz))
990 if (tmp(i:i).eq.' ') then
994 if (nscripts.gt.maxscripts) then
995 call write2log("Number of scripts exceeds maxscripts!")
1000 if (index(trim(tmp),' ').gt.0) then
1001 scripts(i)=tmp(1:index(trim(tmp),' '))
1003 scripts(i)=tmp(1:len_trim(tmp))
1005 call write2log("Script will be used: "//scripts(i))
1006 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
1008 end select ! Koniec czytania "scripts="
1013 c write additional info to log
1014 write(tmp,'(F8.5)') MUTRATIO
1015 call write2log("Mutation ratio is set to "//tmp)
1016 write(tmp,'(I2)') BANK_MULTIPLIER
1017 call write2log("Bank multiplier is set to "//tmp)
1018 write(tmp,'(I5)') BANK_MULTIPLIER*banksize
1020 call write2log("Trail population size will be "//trim(tmp)//" in&
1023 call write2log("Input file unresga.inp does not exist!!")
1024 write(*,*) "Input file unresga.inp does not exist!!"
1028 c Checking set parameters after reading input
1032 call write2log("Can't find 'PDB=' entry in input file")
1036 inquire(FILE=trim(pdbfiles(i)),EXIST=ex)
1038 call write2log("Can't find file "//trim(pdbfiles(i)))
1045 if (maxgen.eq.0) then
1046 call write2log("Couldn't find GENERATIONS= entry in input file")
1051 if (banksize.eq.0) then
1052 call write2log("Can't find POPULATION= entry in input file")
1054 elseif(mod(banksize,2).eq.1) then
1055 call write2log("POPULATION has to be a even number.")
1060 if (ntwham.ne.npdb) then
1061 call write2log("Number of WHAM templates and PDB files is not equ&
1067 if (ntmremd.ne.npdb) then
1068 call write2log("Number of MREMD templates and PDB files is not eq&
1074 if (nscripts.gt.0) then
1076 inquire(FILE=trim(scripts(i)),EXIST=ex)
1078 call write2log("Can't find file "//trim(scripts(i)))
1085 end subroutine ReadInput
1088 c ======================================================================
1089 c CreateInputs subroutine
1090 c ======================================================================
1091 c Creates input files for MREMD, WHAM and ZSCORE.
1092 c ----------------------------------------------------------------------
1094 subroutine CreateInputs(rozmiar,pop)
1097 include 'common.inc'
1098 character(80),dimension(5) :: wagi
1099 integer*4 :: rozmiar,i,j,k
1101 dimension pop(rozmiar,21)
1102 !character*50 :: pdbfiles(maxnpdb)
1103 character*50 :: command,tmptext, prefix
1104 character*256 :: line,plik
1105 integer :: gdzie = 0
1106 character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1107 &8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8&
1109 !write(tmptext,'(I10)') banksize
1110 !call write2log(trim(tmptext))
1113 prefix=pdbfiles(j)(1:len_trim(pdbfiles(j))-4)
1114 c Create MREMD input
1115 open(imremd, file = mremdtemplate(j))
1116 plik=trim(prefix)//"_"//omremdfn
1117 call write2log(trim("Creating MREMD input: "//plik))
1118 open(omremd, file = plik)
1119 3000 read(imremd, '(A)', end = 3010) line
1120 gdzie = index (line, '{PREFIX}')
1121 if (gdzie .ne. 0) then
1122 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1123 write(omremd,"(A)") trim(line)
1124 elseif (index(line,'{POPSIZE}') .ne. 0) then
1125 gdzie = index (line, '{POPSIZE}')
1126 write(tmptext,'(I4)') rozmiar
1127 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(
1129 write(omremd,"(A)") trim(line)
1130 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1131 gdzie = index (line, '{WEIGHTS}')
1133 call Weights2Str(pop(i,:),CUTOFF,wagi)
1134 write(omremd,"(A)") wagi(1)
1135 write(omremd,"(A)") wagi(2)
1136 write(omremd,"(A)") wagi(3)
1137 write(omremd,"(A)") wagi(4)
1138 write(omremd,"(A)") wagi(5)
1141 write(omremd,"(A)") trim(line)
1150 open(iwham, file = whamtemplate(j))
1151 plik=trim(prefix)//"_"//owhamfn
1152 call write2log(trim("Creating WHAM input: "//plik))
1153 open(owham, file = plik)
1154 3015 read(iwham, '(A)', end = 3020) line
1155 gdzie = index (line, '{PREFIX}')
1156 if (gdzie .ne. 0) then
1157 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1158 write(owham,"(A)") trim(line)
1159 elseif (index(line,'{POPSIZE}') .ne. 0) then
1160 gdzie = index (line, '{POPSIZE}')
1161 write(tmptext,'(I4)') rozmiar
1162 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1164 write(owham,"(A)") trim(line)
1165 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1166 gdzie = index (line, '{WEIGHTS}')
1168 call Weights2Str(pop(i,:),CUTOFF,wagi)
1169 write(owham,"(A)") wagi(1)
1170 write(owham,"(A)") wagi(2)
1171 write(owham,"(A)") wagi(3)
1172 write(owham,"(A)") wagi(4)
1173 write(owham,"(A)") wagi(5)
1174 write(owham,"(A)") ''
1177 write(owham,"(A)") trim(line)
1185 c Create ZSCORE input
1186 open(izscore, file = izscorefn)
1187 call write2log(trim("Creating ZSCORE input: "//ozscorefn))
1188 open(ozscore, file = ozscorefn)
1189 3025 read(izscore, '(A)', end = 3030) line
1190 gdzie = index (line, '{PREFIX')
1191 if (gdzie .ne. 0) then
1192 read(line(gdzie+7:gdzie+9),'(I2.2)') k
1193 prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1194 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+10:)&
1197 write(ozscore,"(A)") trim(line)
1198 elseif (index(line,'{PARAM') .ne. 0) then
1199 gdzie=index(line,'{PARAM')
1200 read(line(gdzie+6:gdzie+8),'(I2.2)') k
1201 prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1203 write(tmptext,'(I3.3)') i
1204 line=trim(line(1:gdzie-1))//trim(prefix)//"_wham_par"//tmptex&
1206 write(ozscore,"(A)") trim(line)
1208 elseif (index(line,'{POPSIZE}') .ne. 0) then
1209 gdzie = index (line, '{POPSIZE}')
1210 write(tmptext,'(I4)') rozmiar
1211 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1213 write(ozscore,"(A)") trim(line)
1214 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1215 gdzie = index (line, '{WEIGHTS}')
1217 call Weights2Str(pop(i,:),CUTOFF,wagi)
1218 write(ozscore,"(A)") wagi(1)
1219 write(ozscore,"(A)") wagi(2)
1220 write(ozscore,"(A)") wagi(3)
1221 write(ozscore,"(A)") wagi(4)
1222 write(ozscore,"(A)") trim(wagi(5))
1223 write(ozscore,"(A)") ''
1225 elseif (index(line,'{MAXMIN}') .ne. 0) then
1226 gdzie = index (line, '{MAXMIN}')
1227 write(tmptext,'(I4)') maxminstep
1228 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1230 write(ozscore,"(A)") trim(line)
1232 write(ozscore,"(A)") trim(line)
1241 call write2log('Copying scripts is disabled.')
1243 c ! Open output file
1244 c plik=trim(prefix)//"/"//trim(scripts(i))
1245 c open(oscript, file = plik)
1247 c plik=trim(scripts(i))
1248 c open(iscript, file = plik)
1249 c ! rewrite from input to output
1250 c3030 read(iscript, '(A)', end = 3035) line
1251 c write(oscript,'(A)') trim(line)
1259 c ==========================================
1261 c ==========================================
1264 c write(ow,'(I4)') generation
1265 c write(ow,'(200A)') "# WLONG WSCP WELEC WBOND WANG WSC&
1266 c &LOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC WTURN&
1267 c &3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR"
1269 c write(ow,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),pop&
1270 c &(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),pop&
1271 c &(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,19)
1275 end subroutine CreateInputs
1278 c ======================================================================
1279 c ZnajdzPodobnego function
1280 c ======================================================================
1281 c Finds the most similar individual to osobnik in population(populacja)
1282 c and returns his position in the population
1283 c ----------------------------------------------------------------------
1285 integer*4 function ZnajdzPodobnego(rozmiar,populacja,osobnik,csacu&
1287 integer*4 :: rozmiar,pozycja
1288 real*8 :: delta,mindelta,csacutoff
1290 real*8, dimension(rozmiar,21) :: populacja
1291 real*8, dimension(21) :: osobnik
1293 mindelta = csacutoff
1294 c mindelta = 10000.0
1300 delta=delta+(populacja(i,j)-osobnik(j))**2
1302 c write(tmp,'(2F7.5 )') delta,csacutoff
1303 c call write2log("Delta is "//tmp)
1305 if (delta.lt.mindelta) then
1308 write(tmp,'(F7.5,X,I)') delta,pozycja
1309 call write2log("Delta is "//tmp)
1313 ZnajdzPodobnego = pozycja
1316 c ======================================================================
1317 c WybierzOsobnika function
1318 c ======================================================================
1319 c Selects individual from population using roulette method
1320 c ----------------------------------------------------------------------
1322 integer*4 function WybierzOsobnika(rozmiar,populacja)
1323 integer*4 :: osobnik, rozmiar
1324 real*8 :: los, partsum
1326 real*8, dimension(rozmiar,21) :: populacja
1332 osobnik = osobnik + 1
1333 partsum = partsum + populacja(osobnik,21)
1334 if ((osobnik.eq.rozmiar).or.(partsum.ge.los)) exit
1337 WybierzOsobnika = osobnik
1341 c ======================================================================
1342 c Najlepszy function
1343 c ======================================================================
1344 c Selects the individual with the best score from population
1345 c ----------------------------------------------------------------------
1347 integer*4 function Najlepszy(rozmiar, maks, populacja)
1348 integer*4 :: rozmiar
1350 real*8, dimension(rozmiar,21) :: populacja
1353 if (populacja(i,20).eq.maks) then
1360 c ======================================================================
1361 c WriteState subroutine
1362 c ======================================================================
1363 c Writes current state of calculations in case jobs have to be
1365 c ----------------------------------------------------------------------
1367 subroutine WriteState()
1370 include 'common.inc'
1372 open(ostate,file=ostatefn)
1373 write(ostate,'(I4)') generation
1374 write(ostate,'(L2)') do_optima
1375 write(ostate,'(L2)') do_ga
1376 write(ostate,'(F7.5)') avrd
1378 end subroutine WriteState
1380 c ======================================================================
1381 c ReadPopulation subroutine
1382 c ======================================================================
1383 c Reads population from population.summary file
1384 c ----------------------------------------------------------------------
1386 subroutine ReadPopulation(nind,pop)
1389 character*8 :: dummy
1390 dimension pop(nind,21)
1393 character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1394 &8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8&
1396 read(opopsum,'(A)') dummy
1398 read(opopsum,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),p&
1399 &op(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),p&
1400 &op(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,1&
1403 end subroutine ReadPopulation
1405 c ======================================================================
1406 c ReadOptimaW subroutine
1407 c ======================================================================
1408 c Reads weights from zscore generated weights_opt.zscorepar_XXX files
1410 c nind - number of individuals/files to read
1411 c pop - array where weights are stored
1412 c ----------------------------------------------------------------------
1414 subroutine ReadOptimaW(nind,pop)
1415 integer*4 :: i,nind,stat
1418 character*40 :: filename
1419 character*200 :: tmptxt
1420 dimension pop(nind,21)
1425 c Cosmetic stuff for pretty log
1426 call write2log("==[ Reading weights from zscore files ]=======")
1427 call write2log("# WLONG WSCP WELEC WBOND WANG WS&
1428 &CLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_L&
1429 &OC WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR")
1431 c Actual weight read and error handling
1434 write(tmptxt,'(I3.3)') i
1435 filename = "weights_opt.zscorepar_"//trim(tmptxt)
1436 open(izsoptw, file = filename)
1437 read(izsoptw,"(A4,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)",iosta&
1438 &t=stat) tmptxt,pop(i,1),tmptxt,pop(i,2),tmptxt,pop(i,3),tmptxt,pop&
1439 &(i,4),tmptxt,pop(i,5)
1442 call write2log("-- ERROR while reading "//trim(filename)//"--")
1448 read(izsoptw,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1449 &t=stat) tmptxt,pop(i,6),tmptxt,pop(i,7),tmptxt,pop(i,8),tmptxt,pop&
1450 &(i,9),tmptxt,pop(i,11)
1453 call write2log("-- ERROR while reading "//trim(filename)//"--")
1459 read(izsoptw,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1460 &t=stat) tmptxt,pop(i,12),tmptxt,pop(i,13),tmptxt,pop(i,14),tmptxt,&
1461 &pop(i,15),tmptxt,pop(i,16)
1464 call write2log("-- ERROR while reading "//trim(filename)//"--")
1470 read(izsoptw,"(A7,F7.5,A8,F7.5,A6,F7.5,A)",iostat=stat) tmptxt, &
1471 &pop(i,19),tmptxt,pop(i,17),tmptxt,pop(i,18),tmptxt
1474 call write2log("-- ERROR while reading "//trim(filename)//"--")
1480 read(izsoptw,"(A22,F7.5)",iostat=stat) tmptxt, pop(i,10)
1483 call write2log("-- ERROR while reading "//trim(filename)//"--")
1490 write(tmptxt,'(19F9.5)') (pop(i,j),j=1,19)
1491 call write2log(trim(tmptxt))
1496 end subroutine ReadOptimaW
1499 c ======================================================================
1500 c ReadZEnergy subroutine
1501 c ======================================================================
1502 c Reads Energy from ZSCORE output files zscore_out_parXXX_GB000
1503 c ----------------------------------------------------------------------
1505 subroutine ReadZEnergy(nind,pop)
1509 character*40 :: filename,tmptxt,tmptxt2
1510 character*100 :: wiersz
1511 dimension pop(nind,21)
1514 call write2log("==[ Reading Final Function Value from zscore ]===&
1517 write(tmptxt,'(I3.3)') i
1519 filename=trim(izenergyfn)//trim(tmptxt)//"_GB000"
1520 open(izenergy, file = filename)
1522 read(izenergy, '(A)', iostat=stat) wiersz
1524 if (wiersz(1:22).eq.' Final function value:') then
1525 tmptxt=trim(adjustl(wiersz(23:48)))
1526 read(tmptxt,'(F16.5)') pop(i,20)
1527 write(tmptxt2,'(E15.7)') pop(i,20)
1528 call write2log("Reading score from "//filename//" : "//tmptxt/&
1532 if (pop(i,20).eq.0) then
1533 c Setting FFV to a realy big valiue or otherwise CalcFittness will have problems
1534 pop(i,20)=huge(0.0d0)
1535 call write2log("ERROR while reading FFV from "//filename)
1540 end subroutine ReadZEnergy
1542 c =======================================================================
1543 c WritePopSum subroutine
1544 c =======================================================================
1545 c Writes whole individuals (weights, score, ID, generation) to file
1546 c population.summary This file can be used to plot evolution of the
1547 c weights and scores over time.
1548 c -----------------------------------------------------------------------
1549 subroutine WritePopSum()
1550 logical :: jestplik = .false.
1551 character*5 :: tmptext
1554 c dimension bank(nind,21)
1556 include 'common.inc'
1559 c write(*,*) banksize
1561 write(tmptext,'(I3.3)') generation
1562 inquire(file=opopsumfn, exist=jestplik)
1563 if (.not. jestplik) then
1564 open(opopsum, file = opopsumfn)
1565 write(opopsum, "(A189)") "# WLONG WSCP WELEC WBOND WANG &
1566 & WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC&
1567 & WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR SCO&
1572 open(opopsum, file = opopsumfn, access = "append")
1573 c call ReadPopulation(banksize,populacja)
1574 write(opopsum, "(A189)") "# ==== Generation "//trim(tmptext)//" ==&
1575 &==================================================================&
1576 &==================================================================&
1577 &================================="
1578 c call write2log(banksize)
1581 c write(opopsum,"(F9.5,2F8.5)") populacja(i,1),populacja(i,2),popu&
1583 write(opopsum,"(F9.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,&
1584 &F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F20.5,I8,I8)") bank(i&
1585 &,1),bank(i,2),bank(i,3),bank(i,4),bank(i,5),bank(i,6),bank(i,7),ba&
1586 &nk(i,8),bank(i,9),bank(i,10),bank(i,11),bank(i,12),bank(i,13),bank&
1587 &(i,14),bank(i,15),bank(i,16),bank(i,17),bank(i,18),bank(i,19),bank&
1593 end subroutine WritePopSum
1595 c ======================================================================
1596 c WriteBank subroutine
1597 c ======================================================================
1598 c Writes CSA BANK to file
1599 c ----------------------------------------------------------------------
1600 subroutine WriteBank(b)
1602 include 'common.inc'
1603 real*8,dimension(banksize,21) :: b
1604 character*250 :: tmptext
1605 character*250 :: header
1608 header="# WLONG WSCP WELEC WBOND WANG WSC&
1609 &LOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 &
1610 &WEL_LOC WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCO&
1613 call write2log("Writing Bank to file "//iobankfn)
1614 open(iobank, file = iobankfn)
1615 write(iobank, "(A)") trim(adjustl(header))
1618 write(iobank,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1619 write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1620 call write2log(trim(tmptext))
1625 c ======================================================================
1626 c ReadBank subroutine
1627 c ======================================================================
1628 c Read CSA BANK from file
1629 c ----------------------------------------------------------------------
1630 subroutine ReadBank(b)
1632 include 'common.inc'
1633 real*8,dimension(banksize,21) :: b
1634 character*250 :: tmptext
1636 call write2log("Reading Bank from file "//iobankfn)
1637 open(iobank, file = iobankfn)
1639 read(iobank,"(A)") tmptext
1641 read(iobank,"(19F10.5,E15.7,F10.5)") (b(i,j),j=1,21)
1642 write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1643 call write2log(trim(tmptext))
1648 c ======================================================================
1649 c CalcFitness subroutine
1650 c ======================================================================
1651 c Calculates fitness of every individual in the bank
1652 c ----------------------------------------------------------------------
1653 subroutine CalcFitness(b,fitn)
1654 include 'common.inc'
1655 real*8,dimension(banksize,21) :: b
1656 real*8 :: fitn, sumfitn
1662 fitn=fitn+(1/b(i,20))
1663 sumfitn=sumfitn+b(i,20)
1666 b(i,21)=1/(b(i,20)*fitn)
1667 c b(i,21)=b(i,20)/fitn
1672 c ======================================================================
1673 c CalcAvgDist subroutine
1674 c ======================================================================
1675 c Calculates average distance between individuals in the bank
1676 c ----------------------------------------------------------------------
1677 subroutine CalcAvgDist(b,avgd)
1678 include 'common.inc'
1679 real*8,dimension(banksize,21) :: b
1684 nd = (banksize-1)*banksize/2 ! number of distances to calculate
1689 d=d+(b(i,w)-b(j,w))**2
1696 c -----------------------------------------------------------------------