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 csacutoff=maxco*avrd-generation*avrd*(maxco-minco)/maxgen
205 write(tmptext,'(F7.5)') csacutoff
206 call write2log("CSA cutoff is now set to "//trim(tmptext))
208 do i=1,BANK_MULTIPLIER*banksize
209 write(tmptext,'(I4)') i
210 call write2log("Checking ind "//trim(tmptext))
211 j=ZnajdzPodobnego(banksize,bank,populacja(i,:),csacutoff)
212 if (j.gt.0) then ! W banku jest podobny
213 if (populacja(i,20).lt.bank(j,20)) then
214 write(tmptext,'(I4)') j
215 write(tmptext2,'(I4)') i
216 call write2log("Swaping ind"//trim(tmptext)//" from bank to &
217 &ind "//trim(tmptext2)//" from population")
218 bank(j,:)=populacja(i,:)
220 else ! W banku nie ma podobnego
221 j=FindWorst(banksize,bank)
222 write(tmptext,'(I4)') j
223 write(tmptext2,'(I4)') i
224 if (populacja(i,20).lt.bank(j,20)) then
225 call write2log("Worst in bank is "//trim(tmptext))
226 write(tmptext,'(21F7.5)') bank(j,:)
227 call write2log("BANK:"//trim(tmptext))
228 call write2log("Swaping worst ind in bank to "//trim(tmptext&
230 write(tmptext,'(21F7.5)') populacja(i,:)
231 call write2log("POP :"//trim(tmptext))
232 bank(j,:)=populacja(i,:)
237 c -----------------------------------------------------------------------
238 c Calculate fitness for reproduction
239 c -----------------------------------------------------------------------
240 call CalcFitness(bank,sumfitness)
242 c --- debug begin ---
244 write (*,*) "Dumping bank: "
246 write(*,FMT) bank(i,:)
248 write (*,*) "Sumfitness: ", sumfitness
256 write(*,*) "Some stuff here in the future"
263 c Have we done last generation ?
264 if (generation.eq.maxgen) then
265 call write2log("Maximum number of genarations reached. QUITING.")
268 c Prapering for next generation so increase counter
271 c -----------------------------------------------------------------------
273 c -----------------------------------------------------------------------
274 call write2log("==[ Reproducting ]============================")
277 do i=1,BANK_MULTIPLIER*banksize
278 m=WybierzOsobnika(banksize,bank)
279 write(tmptext,'(I4)') m
280 write(tmptext2,'(I4)') i
281 call write2log("Reproducting individual "//trim(tmptext)//" from&
282 & bank as new individual "//trim(tmptext2))
283 populacja(i,:)=bank(m,:)
285 case('better','betterv2')
286 populacja(1,:)=bank(1,:)
287 call write2log("Reproducting first (best) individual")
288 do i=2,BANK_MULTIPLIER*banksize
289 m=WybierzOsobnika(banksize,bank)
290 write(tmptext,'(I4)') m
291 write(tmptext2,'(I4)') i
292 call write2log("Reproducting individual "//trim(tmptext)//" from&
293 & bank as new individual "//trim(tmptext2))
294 populacja(i,:)=bank(m,:)
298 c -----------------------------------------------------------------------
300 c -----------------------------------------------------------------------
301 do i=1,BANK_MULTIPLIER*banksize
305 do i=1,BANK_MULTIPLIER*banksize
306 if (pairs(i).eq.0) then
307 210 pos=1+int(rand(0)/(1/real(BANK_MULTIPLIER*banksize)))
312 if (pairs(j).eq.pos) then
320 call write2log("==[ Pairing ]==============================")
321 do i=1,BANK_MULTIPLIER*banksize
322 write(tmptext,'(A11,I3.3,A5,I3.3,A21)') "Individual ",i," and ",p&
323 &airs(i)," will be crossed over"
324 call write2log(trim(tmptext))
327 c -----------------------------------------------------------------------
329 c -----------------------------------------------------------------------
330 call write2log("==[ Crossing ]==============================")
331 do i=1,BANK_MULTIPLIER*banksize
332 if (pairs(i) .gt. i) then
333 call CrossoverS(populacja(i,:),populacja(pairs(i),:))
338 call write2log(" WLONG WSCP WELEC WBOND WANG &
339 &WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC WT&
340 &URN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR")
341 do i=1,BANK_MULTIPLIER*banksize
344 write(text,'(A4,I3.3,A2)') "Ind",i,": "
346 write(tmptext,'(F7.5)') populacja(i,j)
347 text = text(1:tmplen) // tmptext
350 call write2log(trim(text))
353 c -----------------------------------------------------------------------
355 c -----------------------------------------------------------------------
357 call write2log("==[ Mutating ]==============================")
359 if (IGNORE_WEIGHT(i).eq.1) then
360 write(tmptext,'(I2)') i
361 call write2log("Parameter "//trim(WNAME(i))//" will be ignored")
364 do i=1,BANK_MULTIPLIER*banksize
365 call MutationV(populacja(i,:),i,MUTRATIO)
367 c for betterv2 and csa restore best individual
369 case('betterv2','csa')
370 populacja(1,:)=bank(1,:)
374 c Genetic Algorithm blokcode stops here
378 c ---------------------------------
379 c Do some final stuff
380 c ---------------------------------
382 c Set the variables before writing to state file
384 if ((.not.do_ga).and.(.not.do_optima)) then
394 do_optima=.not.do_optima
395 if (generation.eq.0) then
408 c Write te current state and population
415 write(tmptext,'(I)') generation+1
416 call write2log("Preparing inputs for next generation ("//trim(adju&
418 call CreateInputs(BANK_MULTIPLIER*banksize,populacja)
420 c All done? Then let's quit.
426 c ======================================================================
428 c MAIN CODE ENDS HERE
430 c ======================================================================
434 2000 deallocate(populacja)
439 c ======================================================================
440 c Clustering subroutines
441 c ======================================================================
443 include 'cluster.inc'
446 c ======================================================================
447 c Write2log subroutine
448 c ======================================================================
450 c ----------------------------------------------------------------------
452 subroutine Write2log(text)
455 character*200 :: tmptext
456 character(len=*) :: text
458 call itime(timeArray)
459 inquire(file=ologfn,exist=ex)
461 open(olog, file=ologfn, access="append")
462 write(olog,"(A)") "==============================================&
464 tmptext = "= UNRES weights optimizer version "//version
465 write(olog,"(A)") tmptext
466 write(olog,"(A)") info
467 write(olog,"(A)") "==============================================&
471 open(olog, file = ologfn , access = "append")
472 write(olog,"(I2.2,A1,I2.2,A1,I2.2,A2,A)") timeArray(1),":",timeAr&
473 &ray(2),":",timeArray(3),": ",text
475 end subroutine Write2log
477 c ======================================================================
478 c GetNBest subroutine
479 c ======================================================================
480 c Fills up the bank population up to banksize with individuals from
481 c inputpop having the best(lowest) score
482 c ----------------------------------------------------------------------
484 subroutine GetNBest(inputpop,bank,banksize)
485 integer*4 :: banksize, i,j,idx
486 real*8 :: inputpop, bank, best, last
488 dimension inputpop(banksize*BANK_MULTIPLIER,21)
489 dimension bank(banksize,21)
495 do i=1,banksize*BANK_MULTIPLIER
496 if ((inputpop(i,20).lt.best).and.(inputpop(i,20).gt.last)) then
501 bank(j,:)=inputpop(idx,:)
504 end subroutine GetNBest
506 c ======================================================================
507 c FindWorst subroutine
508 c ======================================================================
509 c Finds the worst (highest score) individual in population
510 c ----------------------------------------------------------------------
512 integer*4 function FindWorst(rozmiar,pop)
513 integer*4 :: rozmiar, i, idx
515 dimension pop(rozmiar,21)
520 if (pop(i,20).gt.last) then
529 c ======================================================================
530 c GenPopulation subroutine
531 c ======================================================================
532 c Generates a population (pop) of nind individuals with 19 random
533 c wights from a set of ranges (see GA.inc).
534 c ----------------------------------------------------------------------
536 subroutine GenPopulation(nind,pop)
537 integer*4 :: i,j,nind
539 dimension pop(nind,21)
545 pop(j,1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
547 pop(j,2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
549 pop(j,3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
551 pop(j,4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
553 pop(j,5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
555 pop(j,6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
557 pop(j,7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
559 pop(j,8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
561 pop(j,9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
563 pop(j,10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
565 pop(j,11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
567 pop(j,12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
569 pop(j,13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
571 pop(j,14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
573 pop(j,15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
575 pop(j,16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
577 pop(j,17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
579 pop(j,18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
581 pop(j,19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
585 end subroutine GenPopulation
587 c ======================================================================
588 c CrossoverS subroutine
589 c ======================================================================
590 c Does a simple value crossing-over on a pair
591 c of two given individuals, giving two new individuals.
592 c ----------------------------------------------------------------------
594 subroutine CrossoverS(ind1,ind2)
595 real*8 :: ind1(19),ind2(19),temp(19)
597 loc = 1 + int(rand(0)/(1/18.0))
603 end subroutine CrossoverS
605 c ======================================================================
606 c MutationV subroutine
607 c ======================================================================
608 c Does a value mutation for the given individual (ind) with the
609 c given mutation ratio (mratio).
610 c ----------------------------------------------------------------------
612 subroutine MutationV(ind,id,mratio)
618 character*100 :: tmptext
619 character*150 :: logtext
622 if (IGNORE_WEIGHT(i).eq.1) then
623 c write(tmptext,'(I2)') i
624 c call write2log("Skipping weight "//trim(tmptext))
627 if (rand(0)<mratio) then
631 ind(1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
634 ind(2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
637 ind(3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
640 ind(4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
643 ind(5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
646 ind(6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
649 ind(7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
652 ind(8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
655 ind(9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
658 ind(10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
661 ind(11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
664 ind(12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
667 ind(13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
670 ind(14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
673 ind(15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
676 ind(16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
679 ind(17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
682 ind(18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
685 ind(19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
687 write(tmptext,'(I4,A13,I2,A7,F7.5,A5,F7.5)') id," at position "&
688 &,i," from: ",ti," to: ",ind(i)
689 logtext = "Mutation occured in individual "//tmptext
690 call write2log (logtext)
696 c ======================================================================
697 c FunkcjaCeluB function
698 c ======================================================================
699 c Simple (dummy) scoring function
700 c ----------------------------------------------------------------------
702 real*8 function FunkcjaCeluB(osobnik)
704 real*8 :: osobnik(20)
706 FunkcjaCeluB = 20 + -exp(-abs(WLONG-osobnik(1))) &
707 & - exp(-abs(WSCP-osobnik(2))) &
708 & - exp(-abs(WELEC-osobnik(3))) &
709 & - exp(-abs(WBOND-osobnik(4))) &
710 & - exp(-abs(WANG-osobnik(5))) &
711 & - exp(-abs(WSCLOC-osobnik(6))) &
712 & - exp(-abs(WTOR-osobnik(7))) &
713 & - exp(-abs(WTORD-osobnik(8))) &
714 & - exp(-abs(WCORRH-osobnik(9))) &
715 & - exp(-abs(WCORR4-osobnik(10))) &
716 & - exp(-abs(WCORR5-osobnik(11))) &
717 & - exp(-abs(WCORR6-osobnik(12))) &
718 & - exp(-abs(WEL_LOC-osobnik(13))) &
719 & - exp(-abs(WTURN3-osobnik(14))) &
720 & - exp(-abs(WTURN4-osobnik(15))) &
721 & - exp(-abs(WTURN6-osobnik(16))) &
722 & - exp(-abs(WVDWPP-osobnik(17))) &
723 & - exp(-abs(WHPB-osobnik(18))) &
724 & - exp(-abs(WSCCOR-osobnik(19)))
728 c ======================================================================
729 c FunkcjaCeluB2 function
730 c ======================================================================
731 c Simple (dummy) scoring function
732 c ----------------------------------------------------------------------
734 real*8 function FunkcjaCeluB2(osobnik)
736 real*8 :: osobnik(20)
738 FunkcjaCeluB2 = 1 + 4*(WLONG-osobnik(1))**2 &
739 & + 4*(WSCP-osobnik(2))**2 &
740 & + 4*(WELEC-osobnik(3))**2 &
741 & + 4*(WBOND-osobnik(4))**2 &
742 & + 4*(WANG-osobnik(5))**2 &
743 & + 4*(WSCLOC-osobnik(6))**2 &
744 & + 4*(WTOR-osobnik(7))**2 &
745 & + 4*(WTORD-osobnik(8))**2 &
746 & + 4*(WCORRH-osobnik(9))**2 &
747 & + 4*(WCORR4-osobnik(10))**2 &
748 & + 4*(WCORR5-osobnik(11))**2 &
749 & + 4*(WCORR6-osobnik(12))**2 &
750 & + 4*(WEL_LOC-osobnik(13))**2 &
751 & + 4*(WTURN3-osobnik(14))**2 &
752 & + 4*(WTURN4-osobnik(15))**2 &
753 & + 4*(WTURN6-osobnik(16))**2 &
754 & + 4*(WVDWPP-osobnik(17))**2 &
755 & + 4*(WHPB-osobnik(18))**2 &
756 & + 4*(WSCCOR-osobnik(19))**2
760 c ======================================================================
761 c Weights2str subroutine
762 c ======================================================================
763 c Writes weights of individual(osobnik) to string (lista)
764 c ----------------------------------------------------------------------
766 subroutine Weights2Str(osobnik,coff,lista)
767 real*8, dimension(19) :: osobnik
771 character(80) :: lista
779 write(a,"(A6,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)") "WLONG=",os&
780 &obnik(1)," WSCP=",osobnik(2)," WELEC=",osobnik(3)," WBOND=",osobni&
781 &k(4)," WANG=",osobnik(5)," &"
783 write(a,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)") "WSCLOC=",o&
784 &sobnik(6)," WTOR=",osobnik(7)," WTORD=",osobnik(8)," WCORRH=",osob&
785 &nik(9),"WCORR5=",osobnik(11)," &"
787 write(a,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)") "WCORR6=",o&
788 &sobnik(12)," WEL_LOC=",osobnik(13)," WTURN3=",osobnik(14)," WTURN4&
789 &=",osobnik(15),"WTURN6=",osobnik(16)," &"
791 write(a,"(A7,F7.5,A6,F7.5,A8,F7.5,A)") "WVDWPP=",osobnik(17)," WHP&
792 &B=",osobnik(18)," WSCCOR=",osobnik(19)," &
795 write(a,"(A7,F7.5,A8,F7.5)") "CUTOFF=",coff," WCORR4=",osobnik(10)
800 c ======================================================================
801 c FunkcjaCeluR function
802 c ======================================================================
803 c Real scoring function
804 c ----------------------------------------------------------------------
806 real*8 function FunkcjaCeluR(zscore)
808 if (zscore.eq.0) then
811 FunkcjaCeluR = 1/zscore
816 c ======================================================================
817 c ReadInput subroutine
818 c ======================================================================
820 subroutine ReadInput(status)
827 character*100 :: wiersz,tmp
828 inquire(FILE=inpfn,EXIST=ex)
832 call write2log('==[ Reading main input ]======================')
836 read(inp, '(A)', iostat=stat) wiersz
839 if ((wiersz(1:4).eq.'PDB=').or.(wiersz(1:4).eq.'pdb=')) then
841 tmp = wiersz(5:len_trim(wiersz))
843 if (tmp(i:i).eq.' ') then
847 if (npdb.gt.maxnpdb) then
848 call write2log("Number of input PDB exceeds maxnpdb!")
853 if (index(trim(tmp),' ').gt.0) then
854 pdbfiles(i)=tmp(1:index(trim(tmp),' '))
856 pdbfiles(i)=tmp(1:len_trim(tmp))
858 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
860 endif ! Koniec czytania "PDB="
863 if ((wiersz(1:4).eq.'ALG=').or.(wiersz(1:4).eq.'alg=')) then
864 alg=wiersz(5:len_trim(wiersz))
867 call write2log ("Using 'simple' algorithm")
869 call write2log ("Using 'better' algorithm")
871 call write2log ("Using 'betterv2' algorithm")
873 call write2log ("Using 'CSA' algorithm")
875 call write2log ("Using 'Clustering CSA' algorithm")
878 call write2log ("Unknown algorithm. Using 'csa' as default")
883 select case (wiersz(1:12))
884 case('GENERATIONS=','generations=')
885 tmp = wiersz(13:len_trim(wiersz))
886 read(tmp,'(I)') maxgen
890 c select case(wiersz(1:9))
891 c case('CICUTOFF=','cicutoff=')
892 c tmp = wiersz(10:len_trim(wiersz))
893 c read(tmp(1:len_trim(tmp)),'(F7.5)') cicutoff
894 c call write2log("Initial CSA cutoff is set to "//tmp)
898 select case(wiersz(1:6))
899 case('MINCO=','minco=')
900 tmp = wiersz(7:len_trim(wiersz))
901 read(tmp(1:len_trim(tmp)),'(F7.5)') minco
902 call write2log("Minimal CSA cutoff factor is set to "//tmp)
906 select case(wiersz(1:6))
907 case('MAXCO=','maxco=')
908 tmp = wiersz(7:len_trim(wiersz))
909 read(tmp(1:len_trim(tmp)),'(F7.5)') maxco
910 call write2log("Maximal CSA cutoff factor is set to "//tmp)
915 select case(wiersz(1:11))
916 case('POPULATION=','population=')
917 tmp = wiersz(12:len_trim(wiersz))
918 read(tmp,'(I)') banksize
919 call write2log("Bank size is set to "//tmp)
923 select case(wiersz(1:13))
924 case('WHAMTEMPLATE=','whamtemplate=')
925 tmp = wiersz(14:len_trim(wiersz))
928 if (tmp(i:i).eq.' ') then
932 if (ntwham.gt.maxnpdb) then
933 call write2log("Number of WHAM templates exceeds maxnpdb!")
938 if (index(trim(tmp),' ').gt.0) then
939 whamtemplate(i)=tmp(1:index(trim(tmp),' '))
941 whamtemplate(i)=tmp(1:len_trim(tmp))
943 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
948 select case(wiersz(1:14))
949 case('MREMDTEMPLATE=','mremdtemplate=')
950 tmp = wiersz(15:len_trim(wiersz))
953 if (tmp(i:i).eq.' ') then
957 if (ntmremd.gt.maxnpdb) then
958 call write2log("Number of MREMD templates exceeds maxnpdb!")
963 if (index(trim(tmp),' ').gt.0) then
964 mremdtemplate(i)=tmp(1:index(trim(tmp),' '))
966 mremdtemplate(i)=tmp(1:len_trim(tmp))
968 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
973 select case(wiersz(1:7))
974 case('MAXMIN=','maxmin=')
975 tmp = wiersz(8:len_trim(wiersz))
976 read(tmp,'(I)') maxminstep
977 if (maxminstep.gt.0) then
978 call write2log("ZSCORE weights minimalization set to "//tmp)
981 call write2log("ZSCORE weights minimalization disabled")
987 select case(wiersz(1:8))
988 case('SCRIPTS=','scripts=')
990 tmp = wiersz(9:len_trim(wiersz))
992 if (tmp(i:i).eq.' ') then
996 if (nscripts.gt.maxscripts) then
997 call write2log("Number of scripts exceeds maxscripts!")
1002 if (index(trim(tmp),' ').gt.0) then
1003 scripts(i)=tmp(1:index(trim(tmp),' '))
1005 scripts(i)=tmp(1:len_trim(tmp))
1007 call write2log("Script will be used: "//scripts(i))
1008 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
1010 end select ! Koniec czytania "scripts="
1015 c write additional info to log
1016 write(tmp,'(F8.5)') MUTRATIO
1017 call write2log("Mutation ratio is set to "//tmp)
1018 write(tmp,'(I2)') BANK_MULTIPLIER
1019 call write2log("Bank multiplier is set to "//tmp)
1020 write(tmp,'(I5)') BANK_MULTIPLIER*banksize
1022 call write2log("Trail population size will be "//trim(tmp)//" in&
1025 call write2log("Input file unresga.inp does not exist!!")
1026 write(*,*) "Input file unresga.inp does not exist!!"
1030 c Checking set parameters after reading input
1034 call write2log("Can't find 'PDB=' entry in input file")
1038 inquire(FILE=trim(pdbfiles(i)),EXIST=ex)
1040 call write2log("Can't find file "//trim(pdbfiles(i)))
1047 if (maxgen.eq.0) then
1048 call write2log("Couldn't find GENERATIONS= entry in input file")
1053 if (banksize.eq.0) then
1054 call write2log("Can't find POPULATION= entry in input file")
1056 elseif(mod(banksize,2).eq.1) then
1057 call write2log("POPULATION has to be a even number.")
1062 if (ntwham.ne.npdb) then
1063 call write2log("Number of WHAM templates and PDB files is not equ&
1069 if (ntmremd.ne.npdb) then
1070 call write2log("Number of MREMD templates and PDB files is not eq&
1076 if (nscripts.gt.0) then
1078 inquire(FILE=trim(scripts(i)),EXIST=ex)
1080 call write2log("Can't find file "//trim(scripts(i)))
1087 end subroutine ReadInput
1090 c ======================================================================
1091 c CreateInputs subroutine
1092 c ======================================================================
1093 c Creates input files for MREMD, WHAM and ZSCORE.
1094 c ----------------------------------------------------------------------
1096 subroutine CreateInputs(rozmiar,pop)
1099 include 'common.inc'
1100 character(80),dimension(5) :: wagi
1101 integer*4 :: rozmiar,i,j,k
1103 dimension pop(rozmiar,21)
1104 !character*50 :: pdbfiles(maxnpdb)
1105 character*50 :: command,tmptext, prefix
1106 character*256 :: line,plik
1107 integer :: gdzie = 0
1108 character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1109 &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&
1111 !write(tmptext,'(I10)') banksize
1112 !call write2log(trim(tmptext))
1115 prefix=pdbfiles(j)(1:len_trim(pdbfiles(j))-4)
1116 c Create MREMD input
1117 open(imremd, file = mremdtemplate(j))
1118 plik=trim(prefix)//"_"//omremdfn
1119 call write2log(trim("Creating MREMD input: "//plik))
1120 open(omremd, file = plik)
1121 3000 read(imremd, '(A)', end = 3010) line
1122 gdzie = index (line, '{PREFIX}')
1123 if (gdzie .ne. 0) then
1124 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1125 write(omremd,"(A)") trim(line)
1126 elseif (index(line,'{POPSIZE}') .ne. 0) then
1127 gdzie = index (line, '{POPSIZE}')
1128 write(tmptext,'(I4)') rozmiar
1129 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(
1131 write(omremd,"(A)") trim(line)
1132 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1133 gdzie = index (line, '{WEIGHTS}')
1135 call Weights2Str(pop(i,:),CUTOFF,wagi)
1136 write(omremd,"(A)") wagi(1)
1137 write(omremd,"(A)") wagi(2)
1138 write(omremd,"(A)") wagi(3)
1139 write(omremd,"(A)") wagi(4)
1140 write(omremd,"(A)") wagi(5)
1143 write(omremd,"(A)") trim(line)
1152 open(iwham, file = whamtemplate(j))
1153 plik=trim(prefix)//"_"//owhamfn
1154 call write2log(trim("Creating WHAM input: "//plik))
1155 open(owham, file = plik)
1156 3015 read(iwham, '(A)', end = 3020) line
1157 gdzie = index (line, '{PREFIX}')
1158 if (gdzie .ne. 0) then
1159 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1160 write(owham,"(A)") trim(line)
1161 elseif (index(line,'{POPSIZE}') .ne. 0) then
1162 gdzie = index (line, '{POPSIZE}')
1163 write(tmptext,'(I4)') rozmiar
1164 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1166 write(owham,"(A)") trim(line)
1167 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1168 gdzie = index (line, '{WEIGHTS}')
1170 call Weights2Str(pop(i,:),CUTOFF,wagi)
1171 write(owham,"(A)") wagi(1)
1172 write(owham,"(A)") wagi(2)
1173 write(owham,"(A)") wagi(3)
1174 write(owham,"(A)") wagi(4)
1175 write(owham,"(A)") wagi(5)
1176 write(owham,"(A)") ''
1179 write(owham,"(A)") trim(line)
1187 c Create ZSCORE input
1188 open(izscore, file = izscorefn)
1189 call write2log(trim("Creating ZSCORE input: "//ozscorefn))
1190 open(ozscore, file = ozscorefn)
1191 3025 read(izscore, '(A)', end = 3030) line
1192 gdzie = index (line, '{PREFIX')
1193 if (gdzie .ne. 0) then
1194 read(line(gdzie+7:gdzie+9),'(I2.2)') k
1195 prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1196 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+10:)&
1199 write(ozscore,"(A)") trim(line)
1200 elseif (index(line,'{PARAM') .ne. 0) then
1201 gdzie=index(line,'{PARAM')
1202 read(line(gdzie+6:gdzie+8),'(I2.2)') k
1203 prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1205 write(tmptext,'(I3.3)') i
1206 line=trim(line(1:gdzie-1))//trim(prefix)//"_wham_par"//tmptex&
1208 write(ozscore,"(A)") trim(line)
1210 elseif (index(line,'{POPSIZE}') .ne. 0) then
1211 gdzie = index (line, '{POPSIZE}')
1212 write(tmptext,'(I4)') rozmiar
1213 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1215 write(ozscore,"(A)") trim(line)
1216 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1217 gdzie = index (line, '{WEIGHTS}')
1219 call Weights2Str(pop(i,:),CUTOFF,wagi)
1220 write(ozscore,"(A)") wagi(1)
1221 write(ozscore,"(A)") wagi(2)
1222 write(ozscore,"(A)") wagi(3)
1223 write(ozscore,"(A)") wagi(4)
1224 write(ozscore,"(A)") trim(wagi(5))
1225 write(ozscore,"(A)") ''
1227 elseif (index(line,'{MAXMIN}') .ne. 0) then
1228 gdzie = index (line, '{MAXMIN}')
1229 write(tmptext,'(I4)') maxminstep
1230 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1232 write(ozscore,"(A)") trim(line)
1234 write(ozscore,"(A)") trim(line)
1243 call write2log('Copying scripts is disabled.')
1245 c ! Open output file
1246 c plik=trim(prefix)//"/"//trim(scripts(i))
1247 c open(oscript, file = plik)
1249 c plik=trim(scripts(i))
1250 c open(iscript, file = plik)
1251 c ! rewrite from input to output
1252 c3030 read(iscript, '(A)', end = 3035) line
1253 c write(oscript,'(A)') trim(line)
1261 c ==========================================
1263 c ==========================================
1266 c write(ow,'(I4)') generation
1267 c write(ow,'(200A)') "# WLONG WSCP WELEC WBOND WANG WSC&
1268 c &LOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC WTURN&
1269 c &3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR"
1271 c write(ow,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),pop&
1272 c &(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),pop&
1273 c &(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,19)
1277 end subroutine CreateInputs
1280 c ======================================================================
1281 c ZnajdzPodobnego function
1282 c ======================================================================
1283 c Finds the most similar individual to osobnik in population(populacja)
1284 c and returns his position in the population
1285 c ----------------------------------------------------------------------
1287 integer*4 function ZnajdzPodobnego(rozmiar,populacja,osobnik,csacu&
1289 integer*4 :: rozmiar,pozycja
1290 real*8 :: delta,mindelta,csacutoff
1292 real*8, dimension(rozmiar,21) :: populacja
1293 real*8, dimension(21) :: osobnik
1295 mindelta = csacutoff
1296 c mindelta = 10000.0
1302 delta=delta+(populacja(i,j)-osobnik(j))**2
1304 c write(tmp,'(2F7.5 )') delta,csacutoff
1305 c call write2log("Delta is "//tmp)
1307 if (delta.lt.mindelta) then
1310 write(tmp,'(F7.5,X,I)') delta,pozycja
1311 call write2log("Delta is "//tmp)
1315 ZnajdzPodobnego = pozycja
1318 c ======================================================================
1319 c WybierzOsobnika function
1320 c ======================================================================
1321 c Selects individual from population using roulette method
1322 c ----------------------------------------------------------------------
1324 integer*4 function WybierzOsobnika(rozmiar,populacja)
1325 integer*4 :: osobnik, rozmiar
1326 real*8 :: los, partsum
1328 real*8, dimension(rozmiar,21) :: populacja
1334 osobnik = osobnik + 1
1335 partsum = partsum + populacja(osobnik,21)
1336 if ((osobnik.eq.rozmiar).or.(partsum.ge.los)) exit
1339 WybierzOsobnika = osobnik
1343 c ======================================================================
1344 c Najlepszy function
1345 c ======================================================================
1346 c Selects the individual with the best score from population
1347 c ----------------------------------------------------------------------
1349 integer*4 function Najlepszy(rozmiar, maks, populacja)
1350 integer*4 :: rozmiar
1352 real*8, dimension(rozmiar,21) :: populacja
1355 if (populacja(i,20).eq.maks) then
1362 c ======================================================================
1363 c WriteState subroutine
1364 c ======================================================================
1365 c Writes current state of calculations in case jobs have to be
1367 c ----------------------------------------------------------------------
1369 subroutine WriteState()
1372 include 'common.inc'
1374 open(ostate,file=ostatefn)
1375 write(ostate,'(I4)') generation
1376 write(ostate,'(L2)') do_optima
1377 write(ostate,'(L2)') do_ga
1378 write(ostate,'(F7.5)') avrd
1380 end subroutine WriteState
1382 c ======================================================================
1383 c ReadPopulation subroutine
1384 c ======================================================================
1385 c Reads population from population.summary file
1386 c ----------------------------------------------------------------------
1388 subroutine ReadPopulation(nind,pop)
1391 character*8 :: dummy
1392 dimension pop(nind,21)
1395 character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1396 &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&
1398 read(opopsum,'(A)') dummy
1400 read(opopsum,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),p&
1401 &op(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),p&
1402 &op(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,1&
1405 end subroutine ReadPopulation
1407 c ======================================================================
1408 c ReadOptimaW subroutine
1409 c ======================================================================
1410 c Reads weights from zscore generated weights_opt.zscorepar_XXX files
1412 c nind - number of individuals/files to read
1413 c pop - array where weights are stored
1414 c ----------------------------------------------------------------------
1416 subroutine ReadOptimaW(nind,pop)
1417 integer*4 :: i,nind,stat
1420 character*40 :: filename
1421 character*200 :: tmptxt
1422 dimension pop(nind,21)
1427 c Cosmetic stuff for pretty log
1428 call write2log("==[ Reading weights from zscore files ]=======")
1429 call write2log("# WLONG WSCP WELEC WBOND WANG WS&
1430 &CLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_L&
1431 &OC WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR")
1433 c Actual weight read and error handling
1436 write(tmptxt,'(I3.3)') i
1437 filename = "weights_opt.zscorepar_"//trim(tmptxt)
1438 open(izsoptw, file = filename)
1439 read(izsoptw,"(A4,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)",iosta&
1440 &t=stat) tmptxt,pop(i,1),tmptxt,pop(i,2),tmptxt,pop(i,3),tmptxt,pop&
1441 &(i,4),tmptxt,pop(i,5)
1444 call write2log("-- ERROR while reading "//trim(filename)//"--")
1450 read(izsoptw,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1451 &t=stat) tmptxt,pop(i,6),tmptxt,pop(i,7),tmptxt,pop(i,8),tmptxt,pop&
1452 &(i,9),tmptxt,pop(i,11)
1455 call write2log("-- ERROR while reading "//trim(filename)//"--")
1461 read(izsoptw,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1462 &t=stat) tmptxt,pop(i,12),tmptxt,pop(i,13),tmptxt,pop(i,14),tmptxt,&
1463 &pop(i,15),tmptxt,pop(i,16)
1466 call write2log("-- ERROR while reading "//trim(filename)//"--")
1472 read(izsoptw,"(A7,F7.5,A8,F7.5,A6,F7.5,A)",iostat=stat) tmptxt, &
1473 &pop(i,19),tmptxt,pop(i,17),tmptxt,pop(i,18),tmptxt
1476 call write2log("-- ERROR while reading "//trim(filename)//"--")
1482 read(izsoptw,"(A22,F7.5)",iostat=stat) tmptxt, pop(i,10)
1485 call write2log("-- ERROR while reading "//trim(filename)//"--")
1492 write(tmptxt,'(19F9.5)') (pop(i,j),j=1,19)
1493 call write2log(trim(tmptxt))
1498 end subroutine ReadOptimaW
1501 c ======================================================================
1502 c ReadZEnergy subroutine
1503 c ======================================================================
1504 c Reads Energy from ZSCORE output files zscore_out_parXXX_GB000
1505 c ----------------------------------------------------------------------
1507 subroutine ReadZEnergy(nind,pop)
1511 character*40 :: filename,tmptxt,tmptxt2
1512 character*100 :: wiersz
1513 dimension pop(nind,21)
1516 call write2log("==[ Reading Final Function Value from zscore ]===&
1519 write(tmptxt,'(I3.3)') i
1521 filename=trim(izenergyfn)//trim(tmptxt)//"_GB000"
1522 open(izenergy, file = filename)
1524 read(izenergy, '(A)', iostat=stat) wiersz
1526 if (wiersz(1:22).eq.' Final function value:') then
1527 tmptxt=trim(adjustl(wiersz(23:48)))
1528 read(tmptxt,'(F16.5)') pop(i,20)
1529 write(tmptxt2,'(E15.7)') pop(i,20)
1530 call write2log("Reading score from "//filename//" : "//tmptxt/&
1534 if (pop(i,20).eq.0) then
1535 c Setting FFV to a realy big valiue or otherwise CalcFittness will have problems
1536 pop(i,20)=huge(0.0d0)
1537 call write2log("ERROR while reading FFV from "//filename)
1542 end subroutine ReadZEnergy
1544 c =======================================================================
1545 c WritePopSum subroutine
1546 c =======================================================================
1547 c Writes whole individuals (weights, score, ID, generation) to file
1548 c population.summary This file can be used to plot evolution of the
1549 c weights and scores over time.
1550 c -----------------------------------------------------------------------
1551 subroutine WritePopSum()
1552 logical :: jestplik = .false.
1553 character*5 :: tmptext
1556 c dimension bank(nind,21)
1558 include 'common.inc'
1561 c write(*,*) banksize
1563 write(tmptext,'(I3.3)') generation
1564 inquire(file=opopsumfn, exist=jestplik)
1565 if (.not. jestplik) then
1566 open(opopsum, file = opopsumfn)
1567 write(opopsum, "(A189)") "# WLONG WSCP WELEC WBOND WANG &
1568 & WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC&
1569 & WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR SCO&
1574 open(opopsum, file = opopsumfn, access = "append")
1575 c call ReadPopulation(banksize,populacja)
1576 write(opopsum, "(A189)") "# ==== Generation "//trim(tmptext)//" ==&
1577 &==================================================================&
1578 &==================================================================&
1579 &================================="
1580 c call write2log(banksize)
1583 c write(opopsum,"(F9.5,2F8.5)") populacja(i,1),populacja(i,2),popu&
1585 write(opopsum,"(F9.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,&
1586 &F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F20.5,I8,I8)") bank(i&
1587 &,1),bank(i,2),bank(i,3),bank(i,4),bank(i,5),bank(i,6),bank(i,7),ba&
1588 &nk(i,8),bank(i,9),bank(i,10),bank(i,11),bank(i,12),bank(i,13),bank&
1589 &(i,14),bank(i,15),bank(i,16),bank(i,17),bank(i,18),bank(i,19),bank&
1595 end subroutine WritePopSum
1597 c ======================================================================
1598 c WriteBank subroutine
1599 c ======================================================================
1600 c Writes CSA BANK to file
1601 c ----------------------------------------------------------------------
1602 subroutine WriteBank(b)
1604 include 'common.inc'
1605 real*8,dimension(banksize,21) :: b
1606 character*250 :: tmptext
1607 character*250 :: header
1610 header="# WLONG WSCP WELEC WBOND WANG WSC&
1611 &LOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 &
1612 &WEL_LOC WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCO&
1615 call write2log("Writing Bank to file "//iobankfn)
1616 open(iobank, file = iobankfn)
1617 write(iobank, "(A)") trim(adjustl(header))
1620 write(iobank,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1621 write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1622 call write2log(trim(tmptext))
1627 c ======================================================================
1628 c ReadBank subroutine
1629 c ======================================================================
1630 c Read CSA BANK from file
1631 c ----------------------------------------------------------------------
1632 subroutine ReadBank(b)
1634 include 'common.inc'
1635 real*8,dimension(banksize,21) :: b
1636 character*250 :: tmptext
1638 call write2log("Reading Bank from file "//iobankfn)
1639 open(iobank, file = iobankfn)
1641 read(iobank,"(A)") tmptext
1643 read(iobank,"(19F10.5,E15.7,F10.5)") (b(i,j),j=1,21)
1644 write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1645 call write2log(trim(tmptext))
1650 c ======================================================================
1651 c CalcFitness subroutine
1652 c ======================================================================
1653 c Calculates fitness of every individual in the bank
1654 c ----------------------------------------------------------------------
1655 subroutine CalcFitness(b,fitn)
1656 include 'common.inc'
1657 real*8,dimension(banksize,21) :: b
1658 real*8 :: fitn, sumfitn
1664 fitn=fitn+(1/b(i,20))
1665 sumfitn=sumfitn+b(i,20)
1668 b(i,21)=1/(b(i,20)*fitn)
1669 c b(i,21)=b(i,20)/fitn
1674 c ======================================================================
1675 c CalcAvgDist subroutine
1676 c ======================================================================
1677 c Calculates average distance between individuals in the bank
1678 c ----------------------------------------------------------------------
1679 subroutine CalcAvgDist(b,avgd)
1680 include 'common.inc'
1681 real*8,dimension(banksize,21) :: b
1686 nd = (banksize-1)*banksize/2 ! number of distances to calculate
1691 d=d+(b(i,w)-b(j,w))**2
1698 c -----------------------------------------------------------------------