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 c character*500 :: Weights2Str
19 character(80),dimension(5) :: wagi
20 logical :: searching = .false.
21 logical :: logdata = .true.
22 logical :: fileex = .false.
23 logical :: retry = .false.
24 logical :: dir_ex = .false.
25 logical :: debug = .true.
26 character*200 :: text, tmptext,tmptext2
27 character(len=*), parameter :: FMT = "(19F10.5,E15.7,F10.5)"
28 c character(len=*), parameter :: FMT = "(F10.5,F10.5,F10.5,F10.5,F10&
29 c &.5,F10.5,F10.5,F10.5,F10.5,F10.5,F10.5,F10.5,F10.5,F10.5,F10.5,F10&
30 c &.5,F10.5,F10.5,F10.5,F10.5,F10.5)"
32 c =======================================================================
34 c =======================================================================
37 call ReadInput(status)
38 if ((maxminstep.gt.0).and.(generation.eq.1)) then
44 write(*,*) "There ware input file errors - check log."
49 r = rand(timeArray(1)*10000+timeArray(2)*100+timeArray(3))
54 c-------------------------------------------------------
55 c bank & populacja arrays store this
57 c 20 - energy obtained from zscore
59 c-------------------------------------------------------
60 allocate (bank(banksize,21))
61 allocate (populacja(BANK_MULTIPLIER*banksize,21))
62 allocate (pairs(BANK_MULTIPLIER*banksize))
66 do i=1,BANK_MULTIPLIER*banksize
79 inquire(FILE=ostatefn,EXIST=fileex)
80 c No. Prepere next generation
82 open(ostate,file=ostatefn)
83 read(ostate,'(I4)') generation
84 read(ostate,'(L2)') do_optima
85 read(ostate,'(L2)') do_ga
88 call write2log("Doing GA in this step")
90 if (.not.do_optima) then
92 call write2log("ZSCORE weights minimalization disabled for now")
93 generation=generation+1
96 call ReadOptimaW(BANK_MULTIPLIER*banksize,populacja)
98 c Yes. Generate random zero-population
100 call GenPopulation(BANK_MULTIPLIER*banksize,populacja)
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 write(*,*) "cicutoff: ", cicutoff
156 c Fill the bank just after the first time we get the score
158 if (((generation.eq.2).and.(maxminstep.eq.0)).or.((generation.eq&
159 &.2).and.(maxminstep.gt.0))) then
161 c --- debug begin ---
163 write(*,*) "Calling GetNBest..."
166 call GetNBest(populacja,bank,banksize)
167 c --- debug begin ---
170 write(*,FMT) bank(i,:)
175 call CalcFitness(bank,sumfitness)
183 c --- debug begin ---
185 write(*,*) "Szukam podobnego"
191 do i=1,BANK_MULTIPLIER*banksize
192 write(tmptext,'(I4)') i
193 call write2log("Checking ind "//tmptext)
194 j=ZnajdzPodobnego(banksize,bank,populacja(i,:),csacutoff)
196 if (populacja(i,20).lt.bank(j,20)) then
197 write(tmptext,'(I4)') j
198 write(tmptext2,'(I4)') i
199 call write2log("Swaping ind"//trim(tmptext)//" from bank to &
200 &ind "//trim(tmptext2)//" from population")
201 bank(j,:)=populacja(i,:)
204 j=FindWorst(banksize,bank)
205 write(tmptext,'(I4)') j
206 write(tmptext2,'(I4)') i
207 if (populacja(i,20).lt.bank(j,20)) then
208 call write2log("Worst in bank is "//trim(tmptext))
209 call write2log("Swaping worst ind in bank to "//trim(tmptext&
211 bank(j,:)=populacja(i,:)
216 c -----------------------------------------------------------------------
217 c Calculate fitness for reproduction
218 c -----------------------------------------------------------------------
219 call CalcFitness(bank,sumfitness)
221 c --- debug begin ---
223 write (*,*) "Dumping bank: "
225 write(*,FMT) bank(i,:)
227 write (*,*) "Sumfitness: ", sumfitness
232 csacutoff=csacutoff-(cicutoff/maxgen)
233 c csacutoff=cicutoff*(0.8**(iter-1))
235 write(tmptext,'(F7.5)') csacutoff
237 call write2log("CSA cutoff is now set to "//tmptext)
240 write(*,*) "Some stuff here in the future"
247 c Have we done last generation ?
248 if (generation.eq.maxgen) then
249 call write2log("Maximum number of genarations reached. QUITING.")
252 c Prapering for next generation so increase counter
255 c -----------------------------------------------------------------------
257 c -----------------------------------------------------------------------
258 call write2log("==[ Reproducting ]============================")
261 do i=1,BANK_MULTIPLIER*banksize
262 m=WybierzOsobnika(banksize,bank)
263 write(tmptext,'(I4)') m
264 write(tmptext2,'(I4)') i
265 call write2log("Reproducting individual "//trim(tmptext)//" from&
266 & bank as new individual "//trim(tmptext2))
267 populacja(i,:)=bank(m,:)
269 case('better','betterv2')
270 populacja(1,:)=bank(1,:)
271 call write2log("Reproducting first (best) individual")
272 do i=2,BANK_MULTIPLIER*banksize
273 m=WybierzOsobnika(banksize,bank)
274 write(tmptext,'(I4)') m
275 write(tmptext2,'(I4)') i
276 call write2log("Reproducting individual "//trim(tmptext)//" from&
277 & bank as new individual "//trim(tmptext2))
278 populacja(i,:)=bank(m,:)
282 c -----------------------------------------------------------------------
284 c -----------------------------------------------------------------------
285 do i=1,BANK_MULTIPLIER*banksize
289 do i=1,BANK_MULTIPLIER*banksize
290 if (pairs(i).eq.0) then
291 210 pos=1+int(rand(0)/(1/real(BANK_MULTIPLIER*banksize)))
296 if (pairs(j).eq.pos) then
304 call write2log("==[ Pairing ]==============================")
305 do i=1,BANK_MULTIPLIER*banksize
306 write(tmptext,'(A11,I3.3,A5,I3.3,A21)') "Individual ",i," and ",p&
307 &airs(i)," will be crossed over"
308 call write2log(trim(tmptext))
311 c -----------------------------------------------------------------------
313 c -----------------------------------------------------------------------
314 call write2log("==[ Crossing ]==============================")
315 do i=1,BANK_MULTIPLIER*banksize
316 if (pairs(i) .gt. i) then
317 call CrossoverS(populacja(i,:),populacja(pairs(i),:))
322 call write2log(" WLONG WSCP WELEC WBOND WANG &
323 &WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC WT&
324 &URN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR")
325 do i=1,BANK_MULTIPLIER*banksize
328 write(text,'(A4,I3.3,A2)') "Ind",i,": "
330 write(tmptext,'(F7.5)') populacja(i,j)
331 text = text(1:tmplen) // tmptext
334 call write2log(trim(text))
337 c -----------------------------------------------------------------------
339 c -----------------------------------------------------------------------
341 call write2log("==[ Mutating ]==============================")
343 if (IGNORE_WEIGHT(i).eq.1) then
344 write(tmptext,'(I2)') i
345 call write2log("Parameter "//trim(WNAME(i))//" will be ignored")
348 do i=1,BANK_MULTIPLIER*banksize
349 call MutationV(populacja(i,:),i,MUTRATIO)
351 c for betterv2 and csa restore best individual
353 case('betterv2','csa')
354 populacja(1,:)=bank(1,:)
358 c Genetic Algorithm blokcode stops here
362 c ---------------------------------
363 c Do some final stuff
364 c ---------------------------------
366 c Set the variables before writing to state file
368 if ((.not.do_ga).and.(.not.do_optima)) then
378 do_optima=.not.do_optima
379 if (generation.eq.1) then
392 c Write te current state and population
399 write(tmptext,'(I)') generation
400 call write2log("Preparing inputs for generation "//trim(adjustl(tm&
402 call CreateInputs(BANK_MULTIPLIER*banksize,populacja)
404 c All done? Then let's quit.
410 c ======================================================================
412 c MAIN CODE ENDS HERE
414 c ======================================================================
418 2000 deallocate(populacja)
423 c ======================================================================
424 c Clustering subroutines
425 c ======================================================================
427 include 'cluster.inc'
430 c ======================================================================
431 c Write2log subroutine
432 c ======================================================================
434 c ----------------------------------------------------------------------
436 subroutine Write2log(text)
439 character*200 :: tmptext
440 character(len=*) :: text
442 call itime(timeArray)
443 inquire(file=ologfn,exist=ex)
445 open(olog, file=ologfn, access="append")
446 write(olog,"(A)") "==============================================&
448 tmptext = "= UNRES weights optimizer version "//version
449 write(olog,"(A)") tmptext
450 write(olog,"(A)") info
451 write(olog,"(A)") "==============================================&
455 open(olog, file = ologfn , access = "append")
456 write(olog,"(I2.2,A1,I2.2,A1,I2.2,A2,A)") timeArray(1),":",timeAr&
457 &ray(2),":",timeArray(3),": ",text
459 end subroutine Write2log
461 c ======================================================================
462 c GetNBest subroutine
463 c ======================================================================
464 c Fills up the bank population up to banksize with individuals from
465 c inputpop having the best(lowest) score
466 c ----------------------------------------------------------------------
468 subroutine GetNBest(inputpop,bank,banksize)
469 integer*4 :: banksize, i,j,idx
470 real*8 :: inputpop, bank, best, last
472 dimension inputpop(banksize*BANK_MULTIPLIER,21)
473 dimension bank(banksize,21)
478 best=100000000000000000.0
479 do i=1,banksize*BANK_MULTIPLIER
480 if ((inputpop(i,20).lt.best).and.(inputpop(i,20).gt.last)) then
485 bank(j,:)=inputpop(idx,:)
488 end subroutine GetNBest
490 c ======================================================================
491 c FindWorst subroutine
492 c ======================================================================
493 c Finds the worst (highest score) individual in population
494 c ----------------------------------------------------------------------
496 integer*4 function FindWorst(rozmiar,pop)
497 integer*4 :: rozmiar, i, idx
499 dimension pop(rozmiar,21)
504 if (pop(i,20).gt.last) then
513 c ======================================================================
514 c GenPopulation subroutine
515 c ======================================================================
516 c Generates a population (pop) of nind individuals with 19 random
517 c wights from a set of ranges (see GA.inc).
518 c ----------------------------------------------------------------------
520 subroutine GenPopulation(nind,pop)
521 integer*4 :: i,j,nind
523 dimension pop(nind,21)
529 pop(j,1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
531 pop(j,2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
533 pop(j,3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
535 pop(j,4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
537 pop(j,5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
539 pop(j,6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
541 pop(j,7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
543 pop(j,8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
545 pop(j,9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
547 pop(j,10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
549 pop(j,11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
551 pop(j,12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
553 pop(j,13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
555 pop(j,14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
557 pop(j,15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
559 pop(j,16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
561 pop(j,17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
563 pop(j,18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
565 pop(j,19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
569 end subroutine GenPopulation
571 c ======================================================================
572 c CrossoverS subroutine
573 c ======================================================================
574 c Does a simple value crossing-over on a pair
575 c of two given individuals, giving two new individuals.
576 c ----------------------------------------------------------------------
578 subroutine CrossoverS(ind1,ind2)
579 real*8 :: ind1(19),ind2(19),temp(19)
581 loc = 1 + int(rand(0)/(1/18.0))
582 c write (*,*) "Krzyzowanie pomiedzy pozycja ",(loc-1),"a",loc
588 end subroutine CrossoverS
590 c ======================================================================
591 c MutationV subroutine
592 c ======================================================================
593 c Does a value mutation for the given individual (ind) with the
594 c given mutation ratio (mratio).
595 c ----------------------------------------------------------------------
597 subroutine MutationV(ind,id,mratio)
603 character*100 :: tmptext
604 character*150 :: logtext
607 if (IGNORE_WEIGHT(i).eq.1) then
608 c write(tmptext,'(I2)') i
609 c call write2log("Skipping weight "//trim(tmptext))
612 if (rand(0)<mratio) then
616 ind(1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
619 ind(2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
622 ind(3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
625 ind(4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
628 ind(5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
631 ind(6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
634 ind(7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
637 ind(8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
640 ind(9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
643 ind(10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
646 ind(11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
649 ind(12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
652 ind(13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
655 ind(14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
658 ind(15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
661 ind(16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
664 ind(17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
667 ind(18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
670 ind(19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
672 write(tmptext,'(I4,A13,I2,A7,F7.5,A5,F7.5)') id," at position "&
673 &,i," from: ",ti," to: ",ind(i)
674 logtext = "Mutation occured in individual "//tmptext
675 call write2log (logtext)
681 c ======================================================================
682 c FunkcjaCeluB function
683 c ======================================================================
684 c Simple (dummy) scoring function
685 c ----------------------------------------------------------------------
687 real*8 function FunkcjaCeluB(osobnik)
689 real*8 :: osobnik(20)
691 FunkcjaCeluB = 20 + -exp(-abs(WLONG-osobnik(1))) &
692 & - exp(-abs(WSCP-osobnik(2))) &
693 & - exp(-abs(WELEC-osobnik(3))) &
694 & - exp(-abs(WBOND-osobnik(4))) &
695 & - exp(-abs(WANG-osobnik(5))) &
696 & - exp(-abs(WSCLOC-osobnik(6))) &
697 & - exp(-abs(WTOR-osobnik(7))) &
698 & - exp(-abs(WTORD-osobnik(8))) &
699 & - exp(-abs(WCORRH-osobnik(9))) &
700 & - exp(-abs(WCORR4-osobnik(10))) &
701 & - exp(-abs(WCORR5-osobnik(11))) &
702 & - exp(-abs(WCORR6-osobnik(12))) &
703 & - exp(-abs(WEL_LOC-osobnik(13))) &
704 & - exp(-abs(WTURN3-osobnik(14))) &
705 & - exp(-abs(WTURN4-osobnik(15))) &
706 & - exp(-abs(WTURN6-osobnik(16))) &
707 & - exp(-abs(WVDWPP-osobnik(17))) &
708 & - exp(-abs(WHPB-osobnik(18))) &
709 & - exp(-abs(WSCCOR-osobnik(19)))
713 c ======================================================================
714 c FunkcjaCeluB2 function
715 c ======================================================================
716 c Simple (dummy) scoring function
717 c ----------------------------------------------------------------------
719 real*8 function FunkcjaCeluB2(osobnik)
721 real*8 :: osobnik(20)
723 FunkcjaCeluB2 = 1 + 4*(WLONG-osobnik(1))**2 &
724 & + 4*(WSCP-osobnik(2))**2 &
725 & + 4*(WELEC-osobnik(3))**2 &
726 & + 4*(WBOND-osobnik(4))**2 &
727 & + 4*(WANG-osobnik(5))**2 &
728 & + 4*(WSCLOC-osobnik(6))**2 &
729 & + 4*(WTOR-osobnik(7))**2 &
730 & + 4*(WTORD-osobnik(8))**2 &
731 & + 4*(WCORRH-osobnik(9))**2 &
732 & + 4*(WCORR4-osobnik(10))**2 &
733 & + 4*(WCORR5-osobnik(11))**2 &
734 & + 4*(WCORR6-osobnik(12))**2 &
735 & + 4*(WEL_LOC-osobnik(13))**2 &
736 & + 4*(WTURN3-osobnik(14))**2 &
737 & + 4*(WTURN4-osobnik(15))**2 &
738 & + 4*(WTURN6-osobnik(16))**2 &
739 & + 4*(WVDWPP-osobnik(17))**2 &
740 & + 4*(WHPB-osobnik(18))**2 &
741 & + 4*(WSCCOR-osobnik(19))**2
745 c ======================================================================
746 c Weights2str subroutine
747 c ======================================================================
748 c Writes weights of individual(osobnik) to string (lista)
749 c ----------------------------------------------------------------------
751 subroutine Weights2Str(osobnik,coff,lista)
752 real*8, dimension(19) :: osobnik
756 character(80) :: lista
764 write(a,"(A6,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)") "WLONG=",os&
765 &obnik(1)," WSCP=",osobnik(2)," WELEC=",osobnik(3)," WBOND=",osobni&
766 &k(4)," WANG=",osobnik(5)," &"
768 write(a,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)") "WSCLOC=",o&
769 &sobnik(6)," WTOR=",osobnik(7)," WTORD=",osobnik(8)," WCORRH=",osob&
770 &nik(9),"WCORR5=",osobnik(11)," &"
772 write(a,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)") "WCORR6=",o&
773 &sobnik(12)," WEL_LOC=",osobnik(13)," WTURN3=",osobnik(14)," WTURN4&
774 &=",osobnik(15),"WTURN6=",osobnik(16)," &"
776 write(a,"(A7,F7.5,A6,F7.5,A8,F7.5,A)") "WVDWPP=",osobnik(17)," WHP&
777 &B=",osobnik(18)," WSCCOR=",osobnik(19)," &
780 write(a,"(A7,F7.5,A8,F7.5)") "CUTOFF=",coff," WCORR4=",osobnik(10)
785 c ======================================================================
786 c FunkcjaCeluR function
787 c ======================================================================
788 c Real scoring function
789 c ----------------------------------------------------------------------
791 real*8 function FunkcjaCeluR(zscore)
793 if (zscore.eq.0) then
796 FunkcjaCeluR = 1/zscore
801 c ======================================================================
802 c ReadInput subroutine
803 c ======================================================================
805 subroutine ReadInput(status)
812 character*100 :: wiersz,tmp
813 !character*500 :: wiersze = ''
814 inquire(FILE=inpfn,EXIST=ex)
818 call write2log('==[ Reading main input ]======================')
822 read(inp, '(A)', iostat=stat) wiersz
824 if ((wiersz(1:4).eq.'PDB=').or.(wiersz(1:4).eq.'pdb=')) then
826 tmp = wiersz(5:len_trim(wiersz))
828 if (tmp(i:i).eq.' ') then
832 if (npdb.gt.maxnpdb) then
833 call write2log("Number of input PDB exceeds maxnpdb!")
838 if (index(trim(tmp),' ').gt.0) then
839 pdbfiles(i)=tmp(1:index(trim(tmp),' '))
841 pdbfiles(i)=tmp(1:len_trim(tmp))
843 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
845 endif ! Koniec czytania "PDB="
847 if ((wiersz(1:4).eq.'ALG=').or.(wiersz(1:4).eq.'alg=')) then
848 alg=wiersz(5:len_trim(wiersz))
851 call write2log ("Using 'simple' algorithm")
853 call write2log ("Using 'better' algorithm")
855 call write2log ("Using 'betterv2' algorithm")
857 call write2log ("Using 'CSA' algorithm")
859 call write2log ("Using 'Clustering CSA' algorithm")
862 call write2log ("Unknown algorithm. Using 'csa' as default")
864 endif ! Koniec czytania "ALG="
866 select case (wiersz(1:12))
867 case('GENERATIONS=','generations=')
868 tmp = wiersz(13:len_trim(wiersz))
869 read(tmp,'(I)') maxgen
870 end select ! Koniec czytania "GENERATIONS="
873 select case(wiersz(1:9))
874 case('CICUTOFF=','cicutoff=')
875 tmp = wiersz(10:len_trim(wiersz))
876 read(tmp(1:len_trim(tmp)),'(F7.5)') cicutoff
877 call write2log("Initial CSA cutoff is set to "//tmp)
881 select case(wiersz(1:11))
882 case('POPULATION=','population=')
883 tmp = wiersz(12:len_trim(wiersz))
884 read(tmp,'(I)') banksize
885 call write2log("Bank size is set to "//tmp)
889 select case(wiersz(1:13))
890 case('WHAMTEMPLATE=','whamtemplate=')
891 tmp = wiersz(14:len_trim(wiersz))
894 if (tmp(i:i).eq.' ') then
898 if (ntwham.gt.maxnpdb) then
899 call write2log("Number of WHAM templates exceeds maxnpdb!")
904 if (index(trim(tmp),' ').gt.0) then
905 whamtemplate(i)=tmp(1:index(trim(tmp),' '))
907 whamtemplate(i)=tmp(1:len_trim(tmp))
909 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
914 select case(wiersz(1:14))
915 case('MREMDTEMPLATE=','mremdtemplate=')
916 tmp = wiersz(15:len_trim(wiersz))
919 if (tmp(i:i).eq.' ') then
923 if (ntmremd.gt.maxnpdb) then
924 call write2log("Number of MREMD templates exceeds maxnpdb!")
929 if (index(trim(tmp),' ').gt.0) then
930 mremdtemplate(i)=tmp(1:index(trim(tmp),' '))
932 mremdtemplate(i)=tmp(1:len_trim(tmp))
934 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
939 select case(wiersz(1:7))
940 case('MAXMIN=','maxmin=')
941 tmp = wiersz(8:len_trim(wiersz))
942 read(tmp,'(I)') maxminstep
943 if (maxminstep.gt.0) then
944 call write2log("ZSCORE weights minimalization set to "//tmp)
947 call write2log("ZSCORE weights minimalization disabled")
953 select case(wiersz(1:8))
954 case('SCRIPTS=','scripts=')
956 tmp = wiersz(9:len_trim(wiersz))
958 if (tmp(i:i).eq.' ') then
962 if (nscripts.gt.maxscripts) then
963 call write2log("Number of scripts exceeds maxscripts!")
968 if (index(trim(tmp),' ').gt.0) then
969 scripts(i)=tmp(1:index(trim(tmp),' '))
971 scripts(i)=tmp(1:len_trim(tmp))
973 call write2log("Script will be used: "//scripts(i))
974 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
976 end select ! Koniec czytania "scripts="
981 c write additional info to log
982 write(tmp,'(F8.5)') MUTRATIO
983 call write2log("Mutation ratio is set to "//tmp)
984 write(tmp,'(I2)') BANK_MULTIPLIER
985 call write2log("Bank multiplier is set to "//tmp)
986 write(tmp,'(I5)') BANK_MULTIPLIER*banksize
988 call write2log("Trail population size will be "//trim(tmp)//" in&
991 call write2log("Input file unresga.inp does not exist!!")
992 write(*,*) "Input file unresga.inp does not exist!!"
996 c Checking set parameters after reading input
1000 call write2log("Can't find 'PDB=' entry in input file")
1004 inquire(FILE=trim(pdbfiles(i)),EXIST=ex)
1006 call write2log("Can't find file "//trim(pdbfiles(i)))
1013 if (maxgen.eq.0) then
1014 call write2log("Couldn't find GENERATIONS= entry in input file")
1019 if (banksize.eq.0) then
1020 call write2log("Can't find POPULATION= entry in input file")
1022 elseif(mod(banksize,2).eq.1) then
1023 call write2log("POPULATION has to be a even number.")
1028 if (ntwham.ne.npdb) then
1029 call write2log("Number of WHAM templates and PDB files is not equ&
1035 if (ntmremd.ne.npdb) then
1036 call write2log("Number of MREMD templates and PDB files is not eq&
1042 if (nscripts.gt.0) then
1044 inquire(FILE=trim(scripts(i)),EXIST=ex)
1046 call write2log("Can't find file "//trim(scripts(i)))
1053 end subroutine ReadInput
1056 c ======================================================================
1057 c CreateInputs subroutine
1058 c ======================================================================
1059 c Creates input files for MREMD, WHAM and ZSCORE.
1060 c ----------------------------------------------------------------------
1062 subroutine CreateInputs(rozmiar,pop)
1065 include 'common.inc'
1066 character(80),dimension(5) :: wagi
1067 integer*4 :: rozmiar,i,j,k
1069 dimension pop(rozmiar,21)
1070 !character*50 :: pdbfiles(maxnpdb)
1071 character*50 :: command,tmptext, prefix
1072 character*256 :: line,plik
1073 integer :: gdzie = 0
1074 character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1075 &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&
1077 !write(tmptext,'(I10)') banksize
1078 !call write2log(trim(tmptext))
1081 prefix=pdbfiles(j)(1:len_trim(pdbfiles(j))-4)
1082 c Create MREMD input
1083 open(imremd, file = mremdtemplate(j))
1084 plik=trim(prefix)//"_"//omremdfn
1085 call write2log(trim("Creating MREMD input: "//plik))
1086 open(omremd, file = plik)
1087 3000 read(imremd, '(A)', end = 3010) line
1088 gdzie = index (line, '{PREFIX}')
1089 if (gdzie .ne. 0) then
1090 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1091 write(omremd,"(A)") trim(line)
1092 elseif (index(line,'{POPSIZE}') .ne. 0) then
1093 gdzie = index (line, '{POPSIZE}')
1094 write(tmptext,'(I4)') rozmiar
1095 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(
1097 write(omremd,"(A)") trim(line)
1098 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1099 gdzie = index (line, '{WEIGHTS}')
1101 call Weights2Str(pop(i,:),CUTOFF,wagi)
1102 write(omremd,"(A)") wagi(1)
1103 write(omremd,"(A)") wagi(2)
1104 write(omremd,"(A)") wagi(3)
1105 write(omremd,"(A)") wagi(4)
1106 write(omremd,"(A)") wagi(5)
1109 write(omremd,"(A)") trim(line)
1118 open(iwham, file = whamtemplate(j))
1119 plik=trim(prefix)//"_"//owhamfn
1120 call write2log(trim("Creating WHAM input: "//plik))
1121 open(owham, file = plik)
1122 3015 read(iwham, '(A)', end = 3020) line
1123 gdzie = index (line, '{PREFIX}')
1124 if (gdzie .ne. 0) then
1125 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1126 write(owham,"(A)") trim(line)
1127 elseif (index(line,'{POPSIZE}') .ne. 0) then
1128 gdzie = index (line, '{POPSIZE}')
1129 write(tmptext,'(I4)') rozmiar
1130 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1132 write(owham,"(A)") trim(line)
1133 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1134 gdzie = index (line, '{WEIGHTS}')
1136 call Weights2Str(pop(i,:),CUTOFF,wagi)
1137 write(owham,"(A)") wagi(1)
1138 write(owham,"(A)") wagi(2)
1139 write(owham,"(A)") wagi(3)
1140 write(owham,"(A)") wagi(4)
1141 write(owham,"(A)") wagi(5)
1142 write(owham,"(A)") ''
1145 write(owham,"(A)") trim(line)
1153 c Create ZSCORE input
1154 open(izscore, file = izscorefn)
1155 call write2log(trim("Creating ZSCORE input: "//ozscorefn))
1156 open(ozscore, file = ozscorefn)
1157 3025 read(izscore, '(A)', end = 3030) line
1158 gdzie = index (line, '{PREFIX')
1159 if (gdzie .ne. 0) then
1160 read(line(gdzie+7:gdzie+9),'(I2.2)') k
1161 prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1162 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+10:)&
1165 write(ozscore,"(A)") trim(line)
1166 elseif (index(line,'{PARAM') .ne. 0) then
1167 gdzie=index(line,'{PARAM')
1168 read(line(gdzie+6:gdzie+8),'(I2.2)') k
1169 prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1171 write(tmptext,'(I3.3)') i
1172 line=trim(line(1:gdzie-1))//trim(prefix)//"_wham_par"//tmptex&
1174 write(ozscore,"(A)") trim(line)
1176 elseif (index(line,'{POPSIZE}') .ne. 0) then
1177 gdzie = index (line, '{POPSIZE}')
1178 write(tmptext,'(I4)') rozmiar
1179 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1181 write(ozscore,"(A)") trim(line)
1182 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1183 gdzie = index (line, '{WEIGHTS}')
1185 call Weights2Str(pop(i,:),CUTOFF,wagi)
1186 write(ozscore,"(A)") wagi(1)
1187 write(ozscore,"(A)") wagi(2)
1188 write(ozscore,"(A)") wagi(3)
1189 write(ozscore,"(A)") wagi(4)
1190 write(ozscore,"(A)") trim(wagi(5))
1191 write(ozscore,"(A)") ''
1193 elseif (index(line,'{MAXMIN}') .ne. 0) then
1194 gdzie = index (line, '{MAXMIN}')
1195 write(tmptext,'(I4)') maxminstep
1196 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1198 write(ozscore,"(A)") trim(line)
1200 write(ozscore,"(A)") trim(line)
1209 call write2log('Copying scripts is disabled.')
1211 c ! Open output file
1212 c plik=trim(prefix)//"/"//trim(scripts(i))
1213 c open(oscript, file = plik)
1215 c plik=trim(scripts(i))
1216 c open(iscript, file = plik)
1217 c ! rewrite from input to output
1218 c3030 read(iscript, '(A)', end = 3035) line
1219 c write(oscript,'(A)') trim(line)
1227 c ==========================================
1229 c ==========================================
1232 c write(ow,'(I4)') generation
1233 c write(ow,'(200A)') "# WLONG WSCP WELEC WBOND WANG WSC&
1234 c &LOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC WTURN&
1235 c &3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR"
1237 c write(ow,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),pop&
1238 c &(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),pop&
1239 c &(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,19)
1243 c ==========================================
1244 c Testcode before WHAM+ZSCORE integration
1245 c ==========================================
1247 c plik=trim(prefix)//"/"//trim(opopsumfn)
1248 c !print *,trim(plik)
1249 c open(opopsum, file = plik)
1250 c write(opopsum,'(200A)') "# WLONG WSCP WELEC WBOND WANG
1251 c & WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC&
1252 c & WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR SCORE"
1254 c write(opopsum,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5)&
1255 c &,pop(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12)&
1256 c &,pop(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i&
1260 c ==========================================
1263 c command="mkdir zscore"
1264 c call system(command)
1267 end subroutine CreateInputs
1270 c ======================================================================
1271 c ZnajdzPodobnego function
1272 c ======================================================================
1273 c Finds the most similar individual to osobnik in population(populacja)
1274 c and returns his position in the population
1275 c ----------------------------------------------------------------------
1277 integer*4 function ZnajdzPodobnego(rozmiar,populacja,osobnik,csacu&
1279 integer*4 :: rozmiar,pozycja
1280 real*8 :: delta,mindelta,csacutoff
1282 real*8, dimension(rozmiar,21) :: populacja
1283 real*8, dimension(21) :: osobnik
1285 mindelta = csacutoff
1286 c mindelta = 10000.0
1292 delta=delta+(populacja(i,j)-osobnik(j))**2
1294 c write(tmp,'(2F7.5 )') delta,csacutoff
1295 c call write2log("Delta is "//tmp)
1297 if (delta.lt.mindelta) then
1300 write(tmp,'(F7.5,X,I)') delta,pozycja
1301 call write2log("Delta is "//tmp)
1305 ZnajdzPodobnego = pozycja
1308 c ======================================================================
1309 c WybierzOsobnika function
1310 c ======================================================================
1311 c Selects individual from population using roulette method
1312 c ----------------------------------------------------------------------
1314 integer*4 function WybierzOsobnika(rozmiar,populacja)
1315 integer*4 :: osobnik, rozmiar
1316 real*8 :: los, partsum
1318 real*8, dimension(rozmiar,21) :: populacja
1324 osobnik = osobnik + 1
1325 partsum = partsum + populacja(osobnik,21)
1326 if ((osobnik.eq.rozmiar).or.(partsum.ge.los)) exit
1329 WybierzOsobnika = osobnik
1333 c ======================================================================
1334 c Najlepszy function
1335 c ======================================================================
1336 c Selects the individual with the best score from population
1337 c ----------------------------------------------------------------------
1339 integer*4 function Najlepszy(rozmiar, maks, populacja)
1340 integer*4 :: rozmiar
1342 real*8, dimension(rozmiar,21) :: populacja
1345 if (populacja(i,20).eq.maks) then
1352 c ======================================================================
1353 c WriteState subroutine
1354 c ======================================================================
1355 c Writes current state of calculations in case jobs have to be
1357 c ----------------------------------------------------------------------
1359 subroutine WriteState()
1362 include 'common.inc'
1364 open(ostate,file=ostatefn)
1365 write(ostate,'(I4)') generation
1366 write(ostate,'(L2)') do_optima
1367 write(ostate,'(L2)') do_ga
1369 end subroutine WriteState
1371 c ======================================================================
1372 c ReadPopulation subroutine
1373 c ======================================================================
1374 c Reads population from population.summary file
1375 c ----------------------------------------------------------------------
1377 subroutine ReadPopulation(nind,pop)
1380 character*8 :: dummy
1381 dimension pop(nind,21)
1384 character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1385 &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&
1387 read(opopsum,'(A)') dummy
1389 read(opopsum,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),p&
1390 &op(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),p&
1391 &op(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,1&
1394 end subroutine ReadPopulation
1396 c ======================================================================
1397 c ReadOptimaW subroutine
1398 c ======================================================================
1399 c Reads weights from zscore generated weights_opt.zscorepar_XXX files
1401 c nind - number of individuals/files to read
1402 c pop - array where weights are stored
1403 c ----------------------------------------------------------------------
1405 subroutine ReadOptimaW(nind,pop)
1406 integer*4 :: i,nind,stat
1409 character*40 :: filename
1410 character*200 :: tmptxt
1411 dimension pop(nind,21)
1416 c Cosmetic stuff for pretty log
1417 call write2log("==[ Reading weights from zscore files ]=======")
1418 call write2log("# WLONG WSCP WELEC WBOND WANG WS&
1419 &CLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_L&
1420 &OC WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR")
1422 c Actual weight read and error handling
1425 write(tmptxt,'(I3.3)') i
1426 filename = "weights_opt.zscorepar_"//trim(tmptxt)
1427 open(izsoptw, file = filename)
1428 read(izsoptw,"(A4,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)",iosta&
1429 &t=stat) tmptxt,pop(i,1),tmptxt,pop(i,2),tmptxt,pop(i,3),tmptxt,pop&
1430 &(i,4),tmptxt,pop(i,5)
1433 call write2log("-- ERROR while reading "//trim(filename)//"--")
1439 read(izsoptw,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1440 &t=stat) tmptxt,pop(i,6),tmptxt,pop(i,7),tmptxt,pop(i,8),tmptxt,pop&
1441 &(i,9),tmptxt,pop(i,11)
1444 call write2log("-- ERROR while reading "//trim(filename)//"--")
1450 read(izsoptw,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1451 &t=stat) tmptxt,pop(i,12),tmptxt,pop(i,13),tmptxt,pop(i,14),tmptxt,&
1452 &pop(i,15),tmptxt,pop(i,16)
1455 call write2log("-- ERROR while reading "//trim(filename)//"--")
1461 read(izsoptw,"(A7,F7.5,A8,F7.5,A6,F7.5,A)",iostat=stat) tmptxt, &
1462 &pop(i,19),tmptxt,pop(i,17),tmptxt,pop(i,18),tmptxt
1465 call write2log("-- ERROR while reading "//trim(filename)//"--")
1471 read(izsoptw,"(A22,F7.5)",iostat=stat) tmptxt, pop(i,10)
1474 call write2log("-- ERROR while reading "//trim(filename)//"--")
1481 write(tmptxt,'(19F9.5)') (pop(i,j),j=1,19)
1482 call write2log(trim(tmptxt))
1487 end subroutine ReadOptimaW
1490 c ======================================================================
1491 c ReadZEnergy subroutine
1492 c ======================================================================
1493 c Reads Energy from ZSCORE output files zscore_out_parXXX_GB000
1494 c ----------------------------------------------------------------------
1496 subroutine ReadZEnergy(nind,pop)
1500 character*40 :: filename,tmptxt,tmptxt2
1501 character*100 :: wiersz
1502 dimension pop(nind,21)
1505 call write2log("==[ Reading Final Function Value from zscore ]===&
1508 write(tmptxt,'(I3.3)') i
1510 filename=trim(izenergyfn)//trim(tmptxt)//"_GB000"
1511 open(izenergy, file = filename)
1513 read(izenergy, '(A)', iostat=stat) wiersz
1515 if (wiersz(1:22).eq.' Final function value:') then
1516 tmptxt=trim(adjustl(wiersz(23:48)))
1517 read(tmptxt,'(F16.5)') pop(i,20)
1518 write(tmptxt2,'(E15.7)') pop(i,20)
1519 call write2log("Reading score from "//filename//" : "//tmptxt/&
1523 if (pop(i,20).eq.0) then
1524 c Setting FFV to a realy big valiue or otherwise CalcFittness will have problems
1525 pop(i,20)=huge(0.0d0)
1526 call write2log("ERROR while reading FFV from "//filename)
1531 end subroutine ReadZEnergy
1533 c =======================================================================
1534 c WritePopSum subroutine
1535 c =======================================================================
1536 c Writes whole individuals (weights, score, ID, generation) to file
1537 c population.summary This file can be used to plot evolution of the
1538 c weights and scores over time.
1539 c -----------------------------------------------------------------------
1540 subroutine WritePopSum()
1541 logical :: jestplik = .false.
1542 character*5 :: tmptext
1545 c dimension bank(nind,21)
1547 include 'common.inc'
1550 c write(*,*) banksize
1552 write(tmptext,'(I3.3)') generation
1553 inquire(file=opopsumfn, exist=jestplik)
1554 if (.not. jestplik) then
1555 open(opopsum, file = opopsumfn)
1556 write(opopsum, "(A189)") "# WLONG WSCP WELEC WBOND WANG &
1557 & WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC&
1558 & WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR SCO&
1563 open(opopsum, file = opopsumfn, access = "append")
1564 c call ReadPopulation(banksize,populacja)
1565 write(opopsum, "(A189)") "# ==== Generation "//trim(tmptext)//" ==&
1566 &==================================================================&
1567 &==================================================================&
1568 &================================="
1569 c call write2log(banksize)
1572 c write(opopsum,"(F9.5,2F8.5)") populacja(i,1),populacja(i,2),popu&
1574 write(opopsum,"(F9.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,&
1575 &F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F20.5,I8,I8)") bank(i&
1576 &,1),bank(i,2),bank(i,3),bank(i,4),bank(i,5),bank(i,6),bank(i,7),ba&
1577 &nk(i,8),bank(i,9),bank(i,10),bank(i,11),bank(i,12),bank(i,13),bank&
1578 &(i,14),bank(i,15),bank(i,16),bank(i,17),bank(i,18),bank(i,19),bank&
1584 end subroutine WritePopSum
1586 c ======================================================================
1587 c WriteBank subroutine
1588 c ======================================================================
1589 c Writes CSA BANK to file
1590 c ----------------------------------------------------------------------
1591 subroutine WriteBank(b)
1593 include 'common.inc'
1594 real*8,dimension(banksize,21) :: b
1595 character*250 :: tmptext
1596 character*250 :: header
1599 header="# WLONG WSCP WELEC WBOND WANG WSC&
1600 &LOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 &
1601 &WEL_LOC WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCO&
1604 call write2log("Writing Bank to file "//iobankfn)
1605 open(iobank, file = iobankfn)
1606 write(iobank, "(A)") trim(adjustl(header))
1609 write(iobank,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1610 write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1611 call write2log(trim(tmptext))
1616 c ======================================================================
1617 c ReadBank subroutine
1618 c ======================================================================
1619 c Read CSA BANK from file
1620 c ----------------------------------------------------------------------
1621 subroutine ReadBank(b)
1623 include 'common.inc'
1624 real*8,dimension(banksize,21) :: b
1625 character*250 :: tmptext
1627 call write2log("Reading Bank from file "//iobankfn)
1628 open(iobank, file = iobankfn)
1630 read(iobank,"(A)") tmptext
1632 read(iobank,"(19F10.5,E15.7,F10.5)") (b(i,j),j=1,21)
1633 write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1634 call write2log(trim(tmptext))
1639 c ======================================================================
1640 c CalcFitness subroutine
1641 c ======================================================================
1642 c Calculates fitness of every individual in the bank
1643 c ----------------------------------------------------------------------
1644 subroutine CalcFitness(b,fitn)
1645 include 'common.inc'
1646 real*8,dimension(banksize,21) :: b
1647 real*8 :: fitn, sumfitn
1653 fitn=fitn+(1/b(i,20))
1654 sumfitn=sumfitn+b(i,20)
1657 b(i,21)=1/(b(i,20)*fitn)
1658 c b(i,21)=b(i,20)/fitn
1663 c -----------------------------------------------------------------------