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 write(tmptext2,'(21F7.5)') bank(j,:)
219 call write2log("BANK"//trim(tmptext)//":"//trim(tmptext2))
221 write(tmptext,'(I4)') i
222 write(tmptext2,'(21F7.5)') populacja(i,:)
223 call write2log("POP"//trim(tmptext)//":"//trim(tmptext2))
225 bank(j,:)=populacja(i,:)
227 else ! W banku nie ma podobnego
228 j=FindWorst(banksize,bank)
229 write(tmptext,'(I4)') j
230 if (populacja(i,20).lt.bank(j,20)) then
231 call write2log("Worst in bank is "//trim(tmptext))
232 write(tmptext2,'(21F7.5)') bank(j,:)
233 call write2log("BANK"//trim(tmptext)//":"//trim(tmptext2))
235 write(tmptext,'(I4)') i
236 call write2log("Swaping worst ind in bank to "//trim(tmptext&
238 write(tmptext2,'(21F7.5)') populacja(i,:)
239 call write2log("POP"//trim(tmptext)//":"//trim(tmptext2))
240 bank(j,:)=populacja(i,:)
245 c -----------------------------------------------------------------------
246 c Calculate fitness for reproduction
247 c -----------------------------------------------------------------------
248 call CalcFitness(bank,sumfitness)
250 c --- debug begin ---
252 write (*,*) "Dumping bank: "
254 write(*,FMT) bank(i,:)
256 write (*,*) "Sumfitness: ", sumfitness
264 write(*,*) "Some stuff here in the future"
271 c Have we done last generation ?
272 if (generation.eq.maxgen) then
273 call write2log("Maximum number of genarations reached. QUITING.")
276 c Prapering for next generation so increase counter
279 c -----------------------------------------------------------------------
281 c -----------------------------------------------------------------------
282 call write2log("==[ Reproducting ]============================")
285 do i=1,BANK_MULTIPLIER*banksize
286 m=WybierzOsobnika(banksize,bank)
287 write(tmptext,'(I4)') m
288 write(tmptext2,'(I4)') i
289 call write2log("Reproducting individual "//trim(tmptext)//" from&
290 & bank as new individual "//trim(tmptext2))
291 populacja(i,:)=bank(m,:)
293 case('better','betterv2')
294 populacja(1,:)=bank(1,:)
295 call write2log("Reproducting first (best) individual")
296 do i=2,BANK_MULTIPLIER*banksize
297 m=WybierzOsobnika(banksize,bank)
298 write(tmptext,'(I4)') m
299 write(tmptext2,'(I4)') i
300 call write2log("Reproducting individual "//trim(tmptext)//" from&
301 & bank as new individual "//trim(tmptext2))
302 populacja(i,:)=bank(m,:)
306 c -----------------------------------------------------------------------
308 c -----------------------------------------------------------------------
309 do i=1,BANK_MULTIPLIER*banksize
313 do i=1,BANK_MULTIPLIER*banksize
314 if (pairs(i).eq.0) then
315 210 pos=1+int(rand(0)/(1/real(BANK_MULTIPLIER*banksize)))
320 if (pairs(j).eq.pos) then
328 call write2log("==[ Pairing ]==============================")
329 do i=1,BANK_MULTIPLIER*banksize
330 write(tmptext,'(A11,I3.3,A5,I3.3,A21)') "Individual ",i," and ",p&
331 &airs(i)," will be crossed over"
332 call write2log(trim(tmptext))
335 c -----------------------------------------------------------------------
337 c -----------------------------------------------------------------------
338 call write2log("==[ Crossing ]==============================")
339 do i=1,BANK_MULTIPLIER*banksize
340 if (pairs(i) .gt. i) then
341 call CrossoverS(populacja(i,:),populacja(pairs(i),:))
346 call write2log(" WLONG WSCP WELEC WBOND WANG &
347 &WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC WT&
348 &URN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR")
349 do i=1,BANK_MULTIPLIER*banksize
352 write(text,'(A4,I3.3,A2)') "Ind",i,": "
354 write(tmptext,'(F7.5)') populacja(i,j)
355 text = text(1:tmplen) // tmptext
358 call write2log(trim(text))
361 c -----------------------------------------------------------------------
363 c -----------------------------------------------------------------------
365 call write2log("==[ Mutating ]==============================")
367 if (IGNORE_WEIGHT(i).eq.1) then
368 write(tmptext,'(I2)') i
369 call write2log("Parameter "//trim(WNAME(i))//" will be ignored")
372 do i=1,BANK_MULTIPLIER*banksize
373 call MutationV(populacja(i,:),i,MUTRATIO)
375 c for betterv2 and csa restore best individual
377 case('betterv2','csa')
378 populacja(1,:)=bank(1,:)
382 c Genetic Algorithm blokcode stops here
386 c ---------------------------------
387 c Do some final stuff
388 c ---------------------------------
390 c Set the variables before writing to state file
392 if ((.not.do_ga).and.(.not.do_optima)) then
402 do_optima=.not.do_optima
403 if (generation.eq.0) then
416 c Write te current state and population
423 write(tmptext,'(I)') generation+1
424 call write2log("Preparing inputs for next generation ("//trim(adju&
426 call CreateInputs(BANK_MULTIPLIER*banksize,populacja)
428 c All done? Then let's quit.
434 c ======================================================================
436 c MAIN CODE ENDS HERE
438 c ======================================================================
442 2000 deallocate(populacja)
447 c ======================================================================
448 c Clustering subroutines
449 c ======================================================================
451 include 'cluster.inc'
454 c ======================================================================
455 c Write2log subroutine
456 c ======================================================================
458 c ----------------------------------------------------------------------
460 subroutine Write2log(text)
463 character*200 :: tmptext
464 character(len=*) :: text
466 call itime(timeArray)
467 inquire(file=ologfn,exist=ex)
469 open(olog, file=ologfn, access="append")
470 write(olog,"(A)") "==============================================&
472 tmptext = "= UNRES weights optimizer version "//version
473 write(olog,"(A)") tmptext
474 write(olog,"(A)") info
475 write(olog,"(A)") "==============================================&
479 open(olog, file = ologfn , access = "append")
480 write(olog,"(I2.2,A1,I2.2,A1,I2.2,A2,A)") timeArray(1),":",timeAr&
481 &ray(2),":",timeArray(3),": ",text
483 end subroutine Write2log
485 c ======================================================================
486 c GetNBest subroutine
487 c ======================================================================
488 c Fills up the bank population up to banksize with individuals from
489 c inputpop having the best(lowest) score
490 c ----------------------------------------------------------------------
492 subroutine GetNBest(inputpop,bank,banksize)
493 integer*4 :: banksize, i,j,idx
494 real*8 :: inputpop, bank, best, last
496 dimension inputpop(banksize*BANK_MULTIPLIER,21)
497 dimension bank(banksize,21)
503 do i=1,banksize*BANK_MULTIPLIER
504 if ((inputpop(i,20).lt.best).and.(inputpop(i,20).gt.last)) then
509 bank(j,:)=inputpop(idx,:)
512 end subroutine GetNBest
514 c ======================================================================
515 c FindWorst subroutine
516 c ======================================================================
517 c Finds the worst (highest score) individual in population
518 c ----------------------------------------------------------------------
520 integer*4 function FindWorst(rozmiar,pop)
521 integer*4 :: rozmiar, i, idx
523 dimension pop(rozmiar,21)
528 if (pop(i,20).gt.last) then
537 c ======================================================================
538 c GenPopulation subroutine
539 c ======================================================================
540 c Generates a population (pop) of nind individuals with 19 random
541 c wights from a set of ranges (see GA.inc).
542 c ----------------------------------------------------------------------
544 subroutine GenPopulation(nind,pop)
545 integer*4 :: i,j,nind
547 dimension pop(nind,21)
553 pop(j,1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
555 pop(j,2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
557 pop(j,3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
559 pop(j,4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
561 pop(j,5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
563 pop(j,6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
565 pop(j,7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
567 pop(j,8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
569 pop(j,9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
571 pop(j,10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
573 pop(j,11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
575 pop(j,12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
577 pop(j,13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
579 pop(j,14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
581 pop(j,15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
583 pop(j,16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
585 pop(j,17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
587 pop(j,18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
589 pop(j,19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
593 end subroutine GenPopulation
595 c ======================================================================
596 c CrossoverS subroutine
597 c ======================================================================
598 c Does a simple value crossing-over on a pair
599 c of two given individuals, giving two new individuals.
600 c ----------------------------------------------------------------------
602 subroutine CrossoverS(ind1,ind2)
603 real*8 :: ind1(19),ind2(19),temp(19)
605 loc = 1 + int(rand(0)/(1/18.0))
611 end subroutine CrossoverS
613 c ======================================================================
614 c MutationV subroutine
615 c ======================================================================
616 c Does a value mutation for the given individual (ind) with the
617 c given mutation ratio (mratio).
618 c ----------------------------------------------------------------------
620 subroutine MutationV(ind,id,mratio)
626 character*100 :: tmptext
627 character*150 :: logtext
630 if (IGNORE_WEIGHT(i).eq.1) then
631 c write(tmptext,'(I2)') i
632 c call write2log("Skipping weight "//trim(tmptext))
635 if (rand(0)<mratio) then
639 ind(1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
642 ind(2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
645 ind(3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
648 ind(4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
651 ind(5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
654 ind(6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
657 ind(7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
660 ind(8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
663 ind(9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
666 ind(10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
669 ind(11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
672 ind(12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
675 ind(13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
678 ind(14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
681 ind(15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
684 ind(16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
687 ind(17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
690 ind(18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
693 ind(19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
695 write(tmptext,'(I4,A13,I2,A7,F7.5,A5,F7.5)') id," at position "&
696 &,i," from: ",ti," to: ",ind(i)
697 logtext = "Mutation occured in individual "//tmptext
698 call write2log (logtext)
704 c ======================================================================
705 c FunkcjaCeluB function
706 c ======================================================================
707 c Simple (dummy) scoring function
708 c ----------------------------------------------------------------------
710 real*8 function FunkcjaCeluB(osobnik)
712 real*8 :: osobnik(20)
714 FunkcjaCeluB = 20 + -exp(-abs(WLONG-osobnik(1))) &
715 & - exp(-abs(WSCP-osobnik(2))) &
716 & - exp(-abs(WELEC-osobnik(3))) &
717 & - exp(-abs(WBOND-osobnik(4))) &
718 & - exp(-abs(WANG-osobnik(5))) &
719 & - exp(-abs(WSCLOC-osobnik(6))) &
720 & - exp(-abs(WTOR-osobnik(7))) &
721 & - exp(-abs(WTORD-osobnik(8))) &
722 & - exp(-abs(WCORRH-osobnik(9))) &
723 & - exp(-abs(WCORR4-osobnik(10))) &
724 & - exp(-abs(WCORR5-osobnik(11))) &
725 & - exp(-abs(WCORR6-osobnik(12))) &
726 & - exp(-abs(WEL_LOC-osobnik(13))) &
727 & - exp(-abs(WTURN3-osobnik(14))) &
728 & - exp(-abs(WTURN4-osobnik(15))) &
729 & - exp(-abs(WTURN6-osobnik(16))) &
730 & - exp(-abs(WVDWPP-osobnik(17))) &
731 & - exp(-abs(WHPB-osobnik(18))) &
732 & - exp(-abs(WSCCOR-osobnik(19)))
736 c ======================================================================
737 c FunkcjaCeluB2 function
738 c ======================================================================
739 c Simple (dummy) scoring function
740 c ----------------------------------------------------------------------
742 real*8 function FunkcjaCeluB2(osobnik)
744 real*8 :: osobnik(20)
746 FunkcjaCeluB2 = 1 + 4*(WLONG-osobnik(1))**2 &
747 & + 4*(WSCP-osobnik(2))**2 &
748 & + 4*(WELEC-osobnik(3))**2 &
749 & + 4*(WBOND-osobnik(4))**2 &
750 & + 4*(WANG-osobnik(5))**2 &
751 & + 4*(WSCLOC-osobnik(6))**2 &
752 & + 4*(WTOR-osobnik(7))**2 &
753 & + 4*(WTORD-osobnik(8))**2 &
754 & + 4*(WCORRH-osobnik(9))**2 &
755 & + 4*(WCORR4-osobnik(10))**2 &
756 & + 4*(WCORR5-osobnik(11))**2 &
757 & + 4*(WCORR6-osobnik(12))**2 &
758 & + 4*(WEL_LOC-osobnik(13))**2 &
759 & + 4*(WTURN3-osobnik(14))**2 &
760 & + 4*(WTURN4-osobnik(15))**2 &
761 & + 4*(WTURN6-osobnik(16))**2 &
762 & + 4*(WVDWPP-osobnik(17))**2 &
763 & + 4*(WHPB-osobnik(18))**2 &
764 & + 4*(WSCCOR-osobnik(19))**2
768 c ======================================================================
769 c Weights2str subroutine
770 c ======================================================================
771 c Writes weights of individual(osobnik) to string (lista)
772 c ----------------------------------------------------------------------
774 subroutine Weights2Str(osobnik,coff,lista)
775 real*8, dimension(19) :: osobnik
779 character(80) :: lista
787 write(a,"(A6,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)") "WLONG=",os&
788 &obnik(1)," WSCP=",osobnik(2)," WELEC=",osobnik(3)," WBOND=",osobni&
789 &k(4)," WANG=",osobnik(5)," &"
791 write(a,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)") "WSCLOC=",o&
792 &sobnik(6)," WTOR=",osobnik(7)," WTORD=",osobnik(8)," WCORRH=",osob&
793 &nik(9),"WCORR5=",osobnik(11)," &"
795 write(a,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)") "WCORR6=",o&
796 &sobnik(12)," WEL_LOC=",osobnik(13)," WTURN3=",osobnik(14)," WTURN4&
797 &=",osobnik(15),"WTURN6=",osobnik(16)," &"
799 write(a,"(A7,F7.5,A6,F7.5,A8,F7.5,A)") "WVDWPP=",osobnik(17)," WHP&
800 &B=",osobnik(18)," WSCCOR=",osobnik(19)," &
803 write(a,"(A7,F7.5,A8,F7.5)") "CUTOFF=",coff," WCORR4=",osobnik(10)
808 c ======================================================================
809 c FunkcjaCeluR function
810 c ======================================================================
811 c Real scoring function
812 c ----------------------------------------------------------------------
814 real*8 function FunkcjaCeluR(zscore)
816 if (zscore.eq.0) then
819 FunkcjaCeluR = 1/zscore
824 c ======================================================================
825 c ReadInput subroutine
826 c ======================================================================
828 subroutine ReadInput(status)
835 character*100 :: wiersz,tmp
836 inquire(FILE=inpfn,EXIST=ex)
840 call write2log('==[ Reading main input ]======================')
844 read(inp, '(A)', iostat=stat) wiersz
847 if ((wiersz(1:4).eq.'PDB=').or.(wiersz(1:4).eq.'pdb=')) then
849 tmp = wiersz(5:len_trim(wiersz))
851 if (tmp(i:i).eq.' ') then
855 if (npdb.gt.maxnpdb) then
856 call write2log("Number of input PDB exceeds maxnpdb!")
861 if (index(trim(tmp),' ').gt.0) then
862 pdbfiles(i)=tmp(1:index(trim(tmp),' '))
864 pdbfiles(i)=tmp(1:len_trim(tmp))
866 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
868 endif ! Koniec czytania "PDB="
871 if ((wiersz(1:4).eq.'ALG=').or.(wiersz(1:4).eq.'alg=')) then
872 alg=wiersz(5:len_trim(wiersz))
875 call write2log ("Using 'simple' algorithm")
877 call write2log ("Using 'better' algorithm")
879 call write2log ("Using 'betterv2' algorithm")
881 call write2log ("Using 'CSA' algorithm")
883 call write2log ("Using 'Clustering CSA' algorithm")
886 call write2log ("Unknown algorithm. Using 'csa' as default")
891 select case (wiersz(1:12))
892 case('GENERATIONS=','generations=')
893 tmp = wiersz(13:len_trim(wiersz))
894 read(tmp,'(I)') maxgen
898 c select case(wiersz(1:9))
899 c case('CICUTOFF=','cicutoff=')
900 c tmp = wiersz(10:len_trim(wiersz))
901 c read(tmp(1:len_trim(tmp)),'(F7.5)') cicutoff
902 c call write2log("Initial CSA cutoff is set to "//tmp)
906 select case(wiersz(1:6))
907 case('MINCO=','minco=')
908 tmp = wiersz(7:len_trim(wiersz))
909 read(tmp(1:len_trim(tmp)),'(F7.5)') minco
910 call write2log("Minimal CSA cutoff factor is set to "//tmp)
914 select case(wiersz(1:6))
915 case('MAXCO=','maxco=')
916 tmp = wiersz(7:len_trim(wiersz))
917 read(tmp(1:len_trim(tmp)),'(F7.5)') maxco
918 call write2log("Maximal CSA cutoff factor is set to "//tmp)
923 select case(wiersz(1:11))
924 case('POPULATION=','population=')
925 tmp = wiersz(12:len_trim(wiersz))
926 read(tmp,'(I)') banksize
927 call write2log("Bank size is set to "//tmp)
931 select case(wiersz(1:13))
932 case('WHAMTEMPLATE=','whamtemplate=')
933 tmp = wiersz(14:len_trim(wiersz))
936 if (tmp(i:i).eq.' ') then
940 if (ntwham.gt.maxnpdb) then
941 call write2log("Number of WHAM templates exceeds maxnpdb!")
946 if (index(trim(tmp),' ').gt.0) then
947 whamtemplate(i)=tmp(1:index(trim(tmp),' '))
949 whamtemplate(i)=tmp(1:len_trim(tmp))
951 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
956 select case(wiersz(1:14))
957 case('MREMDTEMPLATE=','mremdtemplate=')
958 tmp = wiersz(15:len_trim(wiersz))
961 if (tmp(i:i).eq.' ') then
965 if (ntmremd.gt.maxnpdb) then
966 call write2log("Number of MREMD templates exceeds maxnpdb!")
971 if (index(trim(tmp),' ').gt.0) then
972 mremdtemplate(i)=tmp(1:index(trim(tmp),' '))
974 mremdtemplate(i)=tmp(1:len_trim(tmp))
976 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
981 select case(wiersz(1:7))
982 case('MAXMIN=','maxmin=')
983 tmp = wiersz(8:len_trim(wiersz))
984 read(tmp,'(I)') maxminstep
985 if (maxminstep.gt.0) then
986 call write2log("ZSCORE weights minimalization set to "//tmp)
989 call write2log("ZSCORE weights minimalization disabled")
995 select case(wiersz(1:8))
996 case('SCRIPTS=','scripts=')
998 tmp = wiersz(9:len_trim(wiersz))
1000 if (tmp(i:i).eq.' ') then
1004 if (nscripts.gt.maxscripts) then
1005 call write2log("Number of scripts exceeds maxscripts!")
1010 if (index(trim(tmp),' ').gt.0) then
1011 scripts(i)=tmp(1:index(trim(tmp),' '))
1013 scripts(i)=tmp(1:len_trim(tmp))
1015 call write2log("Script will be used: "//scripts(i))
1016 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
1018 end select ! Koniec czytania "scripts="
1023 c write additional info to log
1024 write(tmp,'(F8.5)') MUTRATIO
1025 call write2log("Mutation ratio is set to "//tmp)
1026 write(tmp,'(I2)') BANK_MULTIPLIER
1027 call write2log("Bank multiplier is set to "//tmp)
1028 write(tmp,'(I5)') BANK_MULTIPLIER*banksize
1030 call write2log("Trail population size will be "//trim(tmp)//" in&
1033 call write2log("Input file unresga.inp does not exist!!")
1034 write(*,*) "Input file unresga.inp does not exist!!"
1038 c Checking set parameters after reading input
1042 call write2log("Can't find 'PDB=' entry in input file")
1046 inquire(FILE=trim(pdbfiles(i)),EXIST=ex)
1048 call write2log("Can't find file "//trim(pdbfiles(i)))
1055 if (maxgen.eq.0) then
1056 call write2log("Couldn't find GENERATIONS= entry in input file")
1061 if (banksize.eq.0) then
1062 call write2log("Can't find POPULATION= entry in input file")
1064 elseif(mod(banksize,2).eq.1) then
1065 call write2log("POPULATION has to be a even number.")
1070 if (ntwham.ne.npdb) then
1071 call write2log("Number of WHAM templates and PDB files is not equ&
1077 if (ntmremd.ne.npdb) then
1078 call write2log("Number of MREMD templates and PDB files is not eq&
1084 if (nscripts.gt.0) then
1086 inquire(FILE=trim(scripts(i)),EXIST=ex)
1088 call write2log("Can't find file "//trim(scripts(i)))
1095 end subroutine ReadInput
1098 c ======================================================================
1099 c CreateInputs subroutine
1100 c ======================================================================
1101 c Creates input files for MREMD, WHAM and ZSCORE.
1102 c ----------------------------------------------------------------------
1104 subroutine CreateInputs(rozmiar,pop)
1107 include 'common.inc'
1108 character(80),dimension(5) :: wagi
1109 integer*4 :: rozmiar,i,j,k
1111 dimension pop(rozmiar,21)
1112 !character*50 :: pdbfiles(maxnpdb)
1113 character*50 :: command,tmptext, prefix
1114 character*256 :: line,plik
1115 integer :: gdzie = 0
1116 character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1117 &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&
1119 !write(tmptext,'(I10)') banksize
1120 !call write2log(trim(tmptext))
1123 prefix=pdbfiles(j)(1:len_trim(pdbfiles(j))-4)
1124 c Create MREMD input
1125 open(imremd, file = mremdtemplate(j))
1126 plik=trim(prefix)//"_"//omremdfn
1127 call write2log(trim("Creating MREMD input: "//plik))
1128 open(omremd, file = plik)
1129 3000 read(imremd, '(A)', end = 3010) line
1130 gdzie = index (line, '{PREFIX}')
1131 if (gdzie .ne. 0) then
1132 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1133 write(omremd,"(A)") trim(line)
1134 elseif (index(line,'{POPSIZE}') .ne. 0) then
1135 gdzie = index (line, '{POPSIZE}')
1136 write(tmptext,'(I4)') rozmiar
1137 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(
1139 write(omremd,"(A)") trim(line)
1140 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1141 gdzie = index (line, '{WEIGHTS}')
1143 call Weights2Str(pop(i,:),CUTOFF,wagi)
1144 write(omremd,"(A)") wagi(1)
1145 write(omremd,"(A)") wagi(2)
1146 write(omremd,"(A)") wagi(3)
1147 write(omremd,"(A)") wagi(4)
1148 write(omremd,"(A)") wagi(5)
1151 write(omremd,"(A)") trim(line)
1160 open(iwham, file = whamtemplate(j))
1161 plik=trim(prefix)//"_"//owhamfn
1162 call write2log(trim("Creating WHAM input: "//plik))
1163 open(owham, file = plik)
1164 3015 read(iwham, '(A)', end = 3020) line
1165 gdzie = index (line, '{PREFIX}')
1166 if (gdzie .ne. 0) then
1167 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1168 write(owham,"(A)") trim(line)
1169 elseif (index(line,'{POPSIZE}') .ne. 0) then
1170 gdzie = index (line, '{POPSIZE}')
1171 write(tmptext,'(I4)') rozmiar
1172 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1174 write(owham,"(A)") trim(line)
1175 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1176 gdzie = index (line, '{WEIGHTS}')
1178 call Weights2Str(pop(i,:),CUTOFF,wagi)
1179 write(owham,"(A)") wagi(1)
1180 write(owham,"(A)") wagi(2)
1181 write(owham,"(A)") wagi(3)
1182 write(owham,"(A)") wagi(4)
1183 write(owham,"(A)") wagi(5)
1184 write(owham,"(A)") ''
1187 write(owham,"(A)") trim(line)
1195 c Create ZSCORE input
1196 open(izscore, file = izscorefn)
1197 call write2log(trim("Creating ZSCORE input: "//ozscorefn))
1198 open(ozscore, file = ozscorefn)
1199 3025 read(izscore, '(A)', end = 3030) line
1200 gdzie = index (line, '{PREFIX')
1201 if (gdzie .ne. 0) then
1202 read(line(gdzie+7:gdzie+9),'(I2.2)') k
1203 prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1204 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+10:)&
1207 write(ozscore,"(A)") trim(line)
1208 elseif (index(line,'{PARAM') .ne. 0) then
1209 gdzie=index(line,'{PARAM')
1210 read(line(gdzie+6:gdzie+8),'(I2.2)') k
1211 prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1213 write(tmptext,'(I3.3)') i
1214 line=trim(line(1:gdzie-1))//trim(prefix)//"_wham_par"//tmptex&
1216 write(ozscore,"(A)") trim(line)
1218 elseif (index(line,'{POPSIZE}') .ne. 0) then
1219 gdzie = index (line, '{POPSIZE}')
1220 write(tmptext,'(I4)') rozmiar
1221 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1223 write(ozscore,"(A)") trim(line)
1224 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1225 gdzie = index (line, '{WEIGHTS}')
1227 call Weights2Str(pop(i,:),CUTOFF,wagi)
1228 write(ozscore,"(A)") wagi(1)
1229 write(ozscore,"(A)") wagi(2)
1230 write(ozscore,"(A)") wagi(3)
1231 write(ozscore,"(A)") wagi(4)
1232 write(ozscore,"(A)") trim(wagi(5))
1233 write(ozscore,"(A)") ''
1235 elseif (index(line,'{MAXMIN}') .ne. 0) then
1236 gdzie = index (line, '{MAXMIN}')
1237 write(tmptext,'(I4)') maxminstep
1238 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1240 write(ozscore,"(A)") trim(line)
1242 write(ozscore,"(A)") trim(line)
1251 call write2log('Copying scripts is disabled.')
1253 c ! Open output file
1254 c plik=trim(prefix)//"/"//trim(scripts(i))
1255 c open(oscript, file = plik)
1257 c plik=trim(scripts(i))
1258 c open(iscript, file = plik)
1259 c ! rewrite from input to output
1260 c3030 read(iscript, '(A)', end = 3035) line
1261 c write(oscript,'(A)') trim(line)
1269 c ==========================================
1271 c ==========================================
1274 c write(ow,'(I4)') generation
1275 c write(ow,'(200A)') "# WLONG WSCP WELEC WBOND WANG WSC&
1276 c &LOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC WTURN&
1277 c &3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR"
1279 c write(ow,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),pop&
1280 c &(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),pop&
1281 c &(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,19)
1285 end subroutine CreateInputs
1288 c ======================================================================
1289 c ZnajdzPodobnego function
1290 c ======================================================================
1291 c Finds the most similar individual to osobnik in population(populacja)
1292 c and returns his position in the population
1293 c ----------------------------------------------------------------------
1295 integer*4 function ZnajdzPodobnego(rozmiar,populacja,osobnik,csacu&
1297 integer*4 :: rozmiar,pozycja
1298 real*8 :: delta,mindelta,csacutoff
1300 real*8, dimension(rozmiar,21) :: populacja
1301 real*8, dimension(21) :: osobnik
1303 mindelta = csacutoff
1304 c mindelta = 10000.0
1310 delta=delta+(populacja(i,j)-osobnik(j))**2
1312 c write(tmp,'(2F7.5 )') delta,csacutoff
1313 c call write2log("Delta is "//tmp)
1315 if (delta.lt.mindelta) then
1318 write(tmp,'(F7.5,X,I)') delta,pozycja
1319 call write2log("Delta is "//tmp)
1323 ZnajdzPodobnego = pozycja
1326 c ======================================================================
1327 c WybierzOsobnika function
1328 c ======================================================================
1329 c Selects individual from population using roulette method
1330 c ----------------------------------------------------------------------
1332 integer*4 function WybierzOsobnika(rozmiar,populacja)
1333 integer*4 :: osobnik, rozmiar
1334 real*8 :: los, partsum
1336 real*8, dimension(rozmiar,21) :: populacja
1342 osobnik = osobnik + 1
1343 partsum = partsum + populacja(osobnik,21)
1344 if ((osobnik.eq.rozmiar).or.(partsum.ge.los)) exit
1347 WybierzOsobnika = osobnik
1351 c ======================================================================
1352 c Najlepszy function
1353 c ======================================================================
1354 c Selects the individual with the best score from population
1355 c ----------------------------------------------------------------------
1357 integer*4 function Najlepszy(rozmiar, maks, populacja)
1358 integer*4 :: rozmiar
1360 real*8, dimension(rozmiar,21) :: populacja
1363 if (populacja(i,20).eq.maks) then
1370 c ======================================================================
1371 c WriteState subroutine
1372 c ======================================================================
1373 c Writes current state of calculations in case jobs have to be
1375 c ----------------------------------------------------------------------
1377 subroutine WriteState()
1380 include 'common.inc'
1382 open(ostate,file=ostatefn)
1383 write(ostate,'(I4)') generation
1384 write(ostate,'(L2)') do_optima
1385 write(ostate,'(L2)') do_ga
1386 write(ostate,'(F7.5)') avrd
1388 end subroutine WriteState
1390 c ======================================================================
1391 c ReadPopulation subroutine
1392 c ======================================================================
1393 c Reads population from population.summary file
1394 c ----------------------------------------------------------------------
1396 subroutine ReadPopulation(nind,pop)
1399 character*8 :: dummy
1400 dimension pop(nind,21)
1403 character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1404 &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&
1406 read(opopsum,'(A)') dummy
1408 read(opopsum,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),p&
1409 &op(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),p&
1410 &op(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,1&
1413 end subroutine ReadPopulation
1415 c ======================================================================
1416 c ReadOptimaW subroutine
1417 c ======================================================================
1418 c Reads weights from zscore generated weights_opt.zscorepar_XXX files
1420 c nind - number of individuals/files to read
1421 c pop - array where weights are stored
1422 c ----------------------------------------------------------------------
1424 subroutine ReadOptimaW(nind,pop)
1425 integer*4 :: i,nind,stat
1428 character*40 :: filename
1429 character*200 :: tmptxt
1430 dimension pop(nind,21)
1435 c Cosmetic stuff for pretty log
1436 call write2log("==[ Reading weights from zscore files ]=======")
1437 call write2log("# WLONG WSCP WELEC WBOND WANG WS&
1438 &CLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_L&
1439 &OC WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR")
1441 c Actual weight read and error handling
1444 write(tmptxt,'(I3.3)') i
1445 filename = "weights_opt.zscorepar_"//trim(tmptxt)
1446 open(izsoptw, file = filename)
1447 read(izsoptw,"(A4,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)",iosta&
1448 &t=stat) tmptxt,pop(i,1),tmptxt,pop(i,2),tmptxt,pop(i,3),tmptxt,pop&
1449 &(i,4),tmptxt,pop(i,5)
1452 call write2log("-- ERROR while reading "//trim(filename)//"--")
1458 read(izsoptw,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1459 &t=stat) tmptxt,pop(i,6),tmptxt,pop(i,7),tmptxt,pop(i,8),tmptxt,pop&
1460 &(i,9),tmptxt,pop(i,11)
1463 call write2log("-- ERROR while reading "//trim(filename)//"--")
1469 read(izsoptw,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1470 &t=stat) tmptxt,pop(i,12),tmptxt,pop(i,13),tmptxt,pop(i,14),tmptxt,&
1471 &pop(i,15),tmptxt,pop(i,16)
1474 call write2log("-- ERROR while reading "//trim(filename)//"--")
1480 read(izsoptw,"(A7,F7.5,A8,F7.5,A6,F7.5,A)",iostat=stat) tmptxt, &
1481 &pop(i,19),tmptxt,pop(i,17),tmptxt,pop(i,18),tmptxt
1484 call write2log("-- ERROR while reading "//trim(filename)//"--")
1490 read(izsoptw,"(A22,F7.5)",iostat=stat) tmptxt, pop(i,10)
1493 call write2log("-- ERROR while reading "//trim(filename)//"--")
1500 write(tmptxt,'(19F9.5)') (pop(i,j),j=1,19)
1501 call write2log(trim(tmptxt))
1506 end subroutine ReadOptimaW
1509 c ======================================================================
1510 c ReadZEnergy subroutine
1511 c ======================================================================
1512 c Reads Energy from ZSCORE output files zscore_out_parXXX_GB000
1513 c ----------------------------------------------------------------------
1515 subroutine ReadZEnergy(nind,pop)
1519 character*40 :: filename,tmptxt,tmptxt2
1520 character*100 :: wiersz
1521 dimension pop(nind,21)
1524 call write2log("==[ Reading Final Function Value from zscore ]===&
1527 write(tmptxt,'(I3.3)') i
1529 filename=trim(izenergyfn)//trim(tmptxt)//"_GB000"
1530 open(izenergy, file = filename)
1532 read(izenergy, '(A)', iostat=stat) wiersz
1534 if (wiersz(1:22).eq.' Final function value:') then
1535 tmptxt=trim(adjustl(wiersz(23:48)))
1536 read(tmptxt,'(F16.5)') pop(i,20)
1537 write(tmptxt2,'(E15.7)') pop(i,20)
1538 call write2log("Reading score from "//filename//" : "//tmptxt/&
1542 if (pop(i,20).eq.0) then
1543 c Setting FFV to a realy big valiue or otherwise CalcFittness will have problems
1544 pop(i,20)=huge(0.0d0)
1545 call write2log("ERROR while reading FFV from "//filename)
1550 end subroutine ReadZEnergy
1552 c =======================================================================
1553 c WritePopSum subroutine
1554 c =======================================================================
1555 c Writes whole individuals (weights, score, ID, generation) to file
1556 c population.summary This file can be used to plot evolution of the
1557 c weights and scores over time.
1558 c -----------------------------------------------------------------------
1559 subroutine WritePopSum()
1560 logical :: jestplik = .false.
1561 character*5 :: tmptext
1564 c dimension bank(nind,21)
1566 include 'common.inc'
1569 c write(*,*) banksize
1571 write(tmptext,'(I3.3)') generation
1572 inquire(file=opopsumfn, exist=jestplik)
1573 if (.not. jestplik) then
1574 open(opopsum, file = opopsumfn)
1575 write(opopsum, "(A189)") "# WLONG WSCP WELEC WBOND WANG &
1576 & WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC&
1577 & WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR SCO&
1582 open(opopsum, file = opopsumfn, access = "append")
1583 c call ReadPopulation(banksize,populacja)
1584 write(opopsum, "(A189)") "# ==== Generation "//trim(tmptext)//" ==&
1585 &==================================================================&
1586 &==================================================================&
1587 &================================="
1588 c call write2log(banksize)
1591 c write(opopsum,"(F9.5,2F8.5)") populacja(i,1),populacja(i,2),popu&
1593 write(opopsum,"(F9.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,&
1594 &F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F20.5,I8,I8)") bank(i&
1595 &,1),bank(i,2),bank(i,3),bank(i,4),bank(i,5),bank(i,6),bank(i,7),ba&
1596 &nk(i,8),bank(i,9),bank(i,10),bank(i,11),bank(i,12),bank(i,13),bank&
1597 &(i,14),bank(i,15),bank(i,16),bank(i,17),bank(i,18),bank(i,19),bank&
1603 end subroutine WritePopSum
1605 c ======================================================================
1606 c WriteBank subroutine
1607 c ======================================================================
1608 c Writes CSA BANK to file
1609 c ----------------------------------------------------------------------
1610 subroutine WriteBank(b)
1612 include 'common.inc'
1613 real*8,dimension(banksize,21) :: b
1614 character*250 :: tmptext
1615 character*250 :: header
1618 header="# WLONG WSCP WELEC WBOND WANG WSC&
1619 &LOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 &
1620 &WEL_LOC WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCO&
1623 call write2log("Writing Bank to file "//iobankfn)
1624 open(iobank, file = iobankfn)
1625 write(iobank, "(A)") trim(adjustl(header))
1628 write(iobank,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1629 write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1630 call write2log(trim(tmptext))
1635 c ======================================================================
1636 c ReadBank subroutine
1637 c ======================================================================
1638 c Read CSA BANK from file
1639 c ----------------------------------------------------------------------
1640 subroutine ReadBank(b)
1642 include 'common.inc'
1643 real*8,dimension(banksize,21) :: b
1644 character*250 :: tmptext
1646 call write2log("Reading Bank from file "//iobankfn)
1647 open(iobank, file = iobankfn)
1649 read(iobank,"(A)") tmptext
1651 read(iobank,"(19F10.5,E15.7,F10.5)") (b(i,j),j=1,21)
1652 write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1653 call write2log(trim(tmptext))
1658 c ======================================================================
1659 c CalcFitness subroutine
1660 c ======================================================================
1661 c Calculates fitness of every individual in the bank
1662 c ----------------------------------------------------------------------
1663 subroutine CalcFitness(b,fitn)
1664 include 'common.inc'
1665 real*8,dimension(banksize,21) :: b
1666 real*8 :: fitn, sumfitn
1672 fitn=fitn+(1/b(i,20))
1673 sumfitn=sumfitn+b(i,20)
1676 b(i,21)=1/(b(i,20)*fitn)
1677 c b(i,21)=b(i,20)/fitn
1682 c ======================================================================
1683 c CalcAvgDist subroutine
1684 c ======================================================================
1685 c Calculates average distance between individuals in the bank
1686 c ----------------------------------------------------------------------
1687 subroutine CalcAvgDist(b,avgd)
1688 include 'common.inc'
1689 real*8,dimension(banksize,21) :: b
1694 nd = (banksize-1)*banksize/2 ! number of distances to calculate
1699 d=d+(b(i,w)-b(j,w))**2
1706 c -----------------------------------------------------------------------