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
84 write(tmptext,'(I4)') generation
85 call write2log("This is genaration "//tmptext)
87 call write2log("Doing GA in this step")
89 if (.not.do_optima) then
91 call write2log("ZSCORE weights minimalization disabled for now")
92 generation=generation+1
95 call ReadOptimaW(BANK_MULTIPLIER*banksize,populacja)
97 c Yes. Generate random zero-population
99 call GenPopulation(BANK_MULTIPLIER*banksize,populacja)
100 write(tmptext,'(I4)') generation
101 call write2log("This is genaration "//tmptext)
106 c End of "fisrst time here?" code
110 c Do we actual use a Genetic Algorithm?
113 if (do_ga.eq..true.) then
115 c Yes. This is only done when doing second pass with ZSCORE (without weight minimalizaton)
116 c or with weight minimalization disabled (maxmin=0) just to obtain final score)
118 c \\---//================================================================
120 c "// Genetic Algorithm code starts here
122 c =//---\\===============================================================
124 call ReadZEnergy(BANK_MULTIPLIER*banksize,populacja)
125 do i=1,BANK_MULTIPLIER*banksize
126 c populacja(i,20)=FunkcjaCeluR(populacja(i,20))
127 write(*,*) "Ind ",i,populacja(i,20)
130 call write2log("==[ Doing GA now ]=============================")
133 c Populate bank depending on algorithm
137 call GetNBest(populacja,bank,banksize)
139 call GetNBest(populacja,bank,banksize)
141 call GetNBest(populacja,bank,banksize)
143 c --- debug begin ---
145 write(*,*) "generation: ",generation
146 write(*,*) "maxminstep: ",maxminstep
147 write(*,*) "do_ga: ",do_ga
148 write(*,*) "do_optima: ",do_optima
149 write(*,*) "do_fs: ",do_fs
150 c write(*,*) "cicutoff: ", cicutoff
157 c Fill the bank just after the first time we get the score
159 if (((generation.eq.1).and.(maxminstep.eq.0)).or.((generation.eq&
160 &.1).and.(maxminstep.gt.0))) then
162 c --- debug begin ---
164 write(*,*) "Calling GetNBest..."
167 call GetNBest(populacja,bank,banksize)
168 call CalcAvgDist(bank,avrd)
169 write(tmptext,'(F7.5)') csacutoff
170 call write2log("CSA cutoff is now set to "//tmptext)
171 csacutoff=(maxco*avrd)-generation*avrd*(maxco-minco)/maxgen
172 c csacutoff=maxco*avrd
176 c --- debug begin ---
179 write(*,FMT) bank(i,:)
184 call CalcFitness(bank,sumfitness)
192 c --- debug begin ---
194 write(*,*) "Szukam podobnego"
199 write(tmptext,'(F7.5)') avrd
200 call write2log("Average distance in bank is "//trim(tmptext))
202 do i=1,BANK_MULTIPLIER*banksize
203 write(tmptext,'(I4)') i
204 call write2log("Checking ind "//tmptext)
205 j=ZnajdzPodobnego(banksize,bank,populacja(i,:),csacutoff)
207 if (populacja(i,20).lt.bank(j,20)) then
208 write(tmptext,'(I4)') j
209 write(tmptext2,'(I4)') i
210 call write2log("Swaping ind"//trim(tmptext)//" from bank to &
211 &ind "//trim(tmptext2)//" from population")
212 bank(j,:)=populacja(i,:)
215 j=FindWorst(banksize,bank)
216 write(tmptext,'(I4)') j
217 write(tmptext2,'(I4)') i
218 if (populacja(i,20).lt.bank(j,20)) then
219 call write2log("Worst in bank is "//trim(tmptext))
220 call write2log("Swaping worst ind in bank to "//trim(tmptext&
222 bank(j,:)=populacja(i,:)
227 c -----------------------------------------------------------------------
228 c Calculate fitness for reproduction
229 c -----------------------------------------------------------------------
230 call CalcFitness(bank,sumfitness)
232 c --- debug begin ---
234 write (*,*) "Dumping bank: "
236 write(*,FMT) bank(i,:)
238 write (*,*) "Sumfitness: ", sumfitness
244 csacutoff=maxco*avrd-generation*avrd*(maxco-minco)/maxgen
245 c csacutoff=csacutoff-(generation*cicutoff/maxgen)
246 c csacutoff=cicutoff*(0.8**(iter-1))
248 write(tmptext,'(F7.5)') csacutoff
249 call write2log("CSA cutoff is now set to "//tmptext)
252 write(*,*) "Some stuff here in the future"
259 c Have we done last generation ?
260 if (generation.eq.maxgen) then
261 call write2log("Maximum number of genarations reached. QUITING.")
264 c Prapering for next generation so increase counter
267 c -----------------------------------------------------------------------
269 c -----------------------------------------------------------------------
270 call write2log("==[ Reproducting ]============================")
273 do i=1,BANK_MULTIPLIER*banksize
274 m=WybierzOsobnika(banksize,bank)
275 write(tmptext,'(I4)') m
276 write(tmptext2,'(I4)') i
277 call write2log("Reproducting individual "//trim(tmptext)//" from&
278 & bank as new individual "//trim(tmptext2))
279 populacja(i,:)=bank(m,:)
281 case('better','betterv2')
282 populacja(1,:)=bank(1,:)
283 call write2log("Reproducting first (best) individual")
284 do i=2,BANK_MULTIPLIER*banksize
285 m=WybierzOsobnika(banksize,bank)
286 write(tmptext,'(I4)') m
287 write(tmptext2,'(I4)') i
288 call write2log("Reproducting individual "//trim(tmptext)//" from&
289 & bank as new individual "//trim(tmptext2))
290 populacja(i,:)=bank(m,:)
294 c -----------------------------------------------------------------------
296 c -----------------------------------------------------------------------
297 do i=1,BANK_MULTIPLIER*banksize
301 do i=1,BANK_MULTIPLIER*banksize
302 if (pairs(i).eq.0) then
303 210 pos=1+int(rand(0)/(1/real(BANK_MULTIPLIER*banksize)))
308 if (pairs(j).eq.pos) then
316 call write2log("==[ Pairing ]==============================")
317 do i=1,BANK_MULTIPLIER*banksize
318 write(tmptext,'(A11,I3.3,A5,I3.3,A21)') "Individual ",i," and ",p&
319 &airs(i)," will be crossed over"
320 call write2log(trim(tmptext))
323 c -----------------------------------------------------------------------
325 c -----------------------------------------------------------------------
326 call write2log("==[ Crossing ]==============================")
327 do i=1,BANK_MULTIPLIER*banksize
328 if (pairs(i) .gt. i) then
329 call CrossoverS(populacja(i,:),populacja(pairs(i),:))
334 call write2log(" WLONG WSCP WELEC WBOND WANG &
335 &WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC WT&
336 &URN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR")
337 do i=1,BANK_MULTIPLIER*banksize
340 write(text,'(A4,I3.3,A2)') "Ind",i,": "
342 write(tmptext,'(F7.5)') populacja(i,j)
343 text = text(1:tmplen) // tmptext
346 call write2log(trim(text))
349 c -----------------------------------------------------------------------
351 c -----------------------------------------------------------------------
353 call write2log("==[ Mutating ]==============================")
355 if (IGNORE_WEIGHT(i).eq.1) then
356 write(tmptext,'(I2)') i
357 call write2log("Parameter "//trim(WNAME(i))//" will be ignored")
360 do i=1,BANK_MULTIPLIER*banksize
361 call MutationV(populacja(i,:),i,MUTRATIO)
363 c for betterv2 and csa restore best individual
365 case('betterv2','csa')
366 populacja(1,:)=bank(1,:)
370 c Genetic Algorithm blokcode stops here
374 c ---------------------------------
375 c Do some final stuff
376 c ---------------------------------
378 c Set the variables before writing to state file
380 if ((.not.do_ga).and.(.not.do_optima)) then
390 do_optima=.not.do_optima
391 if (generation.eq.1) then
404 c Write te current state and population
411 write(tmptext,'(I)') generation
412 call write2log("Preparing inputs for next generation ("//trim(adju&
414 call CreateInputs(BANK_MULTIPLIER*banksize,populacja)
416 c All done? Then let's quit.
422 c ======================================================================
424 c MAIN CODE ENDS HERE
426 c ======================================================================
430 2000 deallocate(populacja)
435 c ======================================================================
436 c Clustering subroutines
437 c ======================================================================
439 include 'cluster.inc'
442 c ======================================================================
443 c Write2log subroutine
444 c ======================================================================
446 c ----------------------------------------------------------------------
448 subroutine Write2log(text)
451 character*200 :: tmptext
452 character(len=*) :: text
454 call itime(timeArray)
455 inquire(file=ologfn,exist=ex)
457 open(olog, file=ologfn, access="append")
458 write(olog,"(A)") "==============================================&
460 tmptext = "= UNRES weights optimizer version "//version
461 write(olog,"(A)") tmptext
462 write(olog,"(A)") info
463 write(olog,"(A)") "==============================================&
467 open(olog, file = ologfn , access = "append")
468 write(olog,"(I2.2,A1,I2.2,A1,I2.2,A2,A)") timeArray(1),":",timeAr&
469 &ray(2),":",timeArray(3),": ",text
471 end subroutine Write2log
473 c ======================================================================
474 c GetNBest subroutine
475 c ======================================================================
476 c Fills up the bank population up to banksize with individuals from
477 c inputpop having the best(lowest) score
478 c ----------------------------------------------------------------------
480 subroutine GetNBest(inputpop,bank,banksize)
481 integer*4 :: banksize, i,j,idx
482 real*8 :: inputpop, bank, best, last
484 dimension inputpop(banksize*BANK_MULTIPLIER,21)
485 dimension bank(banksize,21)
491 do i=1,banksize*BANK_MULTIPLIER
492 if ((inputpop(i,20).lt.best).and.(inputpop(i,20).gt.last)) then
497 bank(j,:)=inputpop(idx,:)
500 end subroutine GetNBest
502 c ======================================================================
503 c FindWorst subroutine
504 c ======================================================================
505 c Finds the worst (highest score) individual in population
506 c ----------------------------------------------------------------------
508 integer*4 function FindWorst(rozmiar,pop)
509 integer*4 :: rozmiar, i, idx
511 dimension pop(rozmiar,21)
516 if (pop(i,20).gt.last) then
525 c ======================================================================
526 c GenPopulation subroutine
527 c ======================================================================
528 c Generates a population (pop) of nind individuals with 19 random
529 c wights from a set of ranges (see GA.inc).
530 c ----------------------------------------------------------------------
532 subroutine GenPopulation(nind,pop)
533 integer*4 :: i,j,nind
535 dimension pop(nind,21)
541 pop(j,1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
543 pop(j,2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
545 pop(j,3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
547 pop(j,4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
549 pop(j,5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
551 pop(j,6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
553 pop(j,7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
555 pop(j,8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
557 pop(j,9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
559 pop(j,10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
561 pop(j,11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
563 pop(j,12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
565 pop(j,13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
567 pop(j,14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
569 pop(j,15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
571 pop(j,16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
573 pop(j,17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
575 pop(j,18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
577 pop(j,19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
581 end subroutine GenPopulation
583 c ======================================================================
584 c CrossoverS subroutine
585 c ======================================================================
586 c Does a simple value crossing-over on a pair
587 c of two given individuals, giving two new individuals.
588 c ----------------------------------------------------------------------
590 subroutine CrossoverS(ind1,ind2)
591 real*8 :: ind1(19),ind2(19),temp(19)
593 loc = 1 + int(rand(0)/(1/18.0))
599 end subroutine CrossoverS
601 c ======================================================================
602 c MutationV subroutine
603 c ======================================================================
604 c Does a value mutation for the given individual (ind) with the
605 c given mutation ratio (mratio).
606 c ----------------------------------------------------------------------
608 subroutine MutationV(ind,id,mratio)
614 character*100 :: tmptext
615 character*150 :: logtext
618 if (IGNORE_WEIGHT(i).eq.1) then
619 c write(tmptext,'(I2)') i
620 c call write2log("Skipping weight "//trim(tmptext))
623 if (rand(0)<mratio) then
627 ind(1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
630 ind(2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
633 ind(3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
636 ind(4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
639 ind(5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
642 ind(6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
645 ind(7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
648 ind(8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
651 ind(9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
654 ind(10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
657 ind(11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
660 ind(12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
663 ind(13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
666 ind(14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
669 ind(15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
672 ind(16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
675 ind(17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
678 ind(18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
681 ind(19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
683 write(tmptext,'(I4,A13,I2,A7,F7.5,A5,F7.5)') id," at position "&
684 &,i," from: ",ti," to: ",ind(i)
685 logtext = "Mutation occured in individual "//tmptext
686 call write2log (logtext)
692 c ======================================================================
693 c FunkcjaCeluB function
694 c ======================================================================
695 c Simple (dummy) scoring function
696 c ----------------------------------------------------------------------
698 real*8 function FunkcjaCeluB(osobnik)
700 real*8 :: osobnik(20)
702 FunkcjaCeluB = 20 + -exp(-abs(WLONG-osobnik(1))) &
703 & - exp(-abs(WSCP-osobnik(2))) &
704 & - exp(-abs(WELEC-osobnik(3))) &
705 & - exp(-abs(WBOND-osobnik(4))) &
706 & - exp(-abs(WANG-osobnik(5))) &
707 & - exp(-abs(WSCLOC-osobnik(6))) &
708 & - exp(-abs(WTOR-osobnik(7))) &
709 & - exp(-abs(WTORD-osobnik(8))) &
710 & - exp(-abs(WCORRH-osobnik(9))) &
711 & - exp(-abs(WCORR4-osobnik(10))) &
712 & - exp(-abs(WCORR5-osobnik(11))) &
713 & - exp(-abs(WCORR6-osobnik(12))) &
714 & - exp(-abs(WEL_LOC-osobnik(13))) &
715 & - exp(-abs(WTURN3-osobnik(14))) &
716 & - exp(-abs(WTURN4-osobnik(15))) &
717 & - exp(-abs(WTURN6-osobnik(16))) &
718 & - exp(-abs(WVDWPP-osobnik(17))) &
719 & - exp(-abs(WHPB-osobnik(18))) &
720 & - exp(-abs(WSCCOR-osobnik(19)))
724 c ======================================================================
725 c FunkcjaCeluB2 function
726 c ======================================================================
727 c Simple (dummy) scoring function
728 c ----------------------------------------------------------------------
730 real*8 function FunkcjaCeluB2(osobnik)
732 real*8 :: osobnik(20)
734 FunkcjaCeluB2 = 1 + 4*(WLONG-osobnik(1))**2 &
735 & + 4*(WSCP-osobnik(2))**2 &
736 & + 4*(WELEC-osobnik(3))**2 &
737 & + 4*(WBOND-osobnik(4))**2 &
738 & + 4*(WANG-osobnik(5))**2 &
739 & + 4*(WSCLOC-osobnik(6))**2 &
740 & + 4*(WTOR-osobnik(7))**2 &
741 & + 4*(WTORD-osobnik(8))**2 &
742 & + 4*(WCORRH-osobnik(9))**2 &
743 & + 4*(WCORR4-osobnik(10))**2 &
744 & + 4*(WCORR5-osobnik(11))**2 &
745 & + 4*(WCORR6-osobnik(12))**2 &
746 & + 4*(WEL_LOC-osobnik(13))**2 &
747 & + 4*(WTURN3-osobnik(14))**2 &
748 & + 4*(WTURN4-osobnik(15))**2 &
749 & + 4*(WTURN6-osobnik(16))**2 &
750 & + 4*(WVDWPP-osobnik(17))**2 &
751 & + 4*(WHPB-osobnik(18))**2 &
752 & + 4*(WSCCOR-osobnik(19))**2
756 c ======================================================================
757 c Weights2str subroutine
758 c ======================================================================
759 c Writes weights of individual(osobnik) to string (lista)
760 c ----------------------------------------------------------------------
762 subroutine Weights2Str(osobnik,coff,lista)
763 real*8, dimension(19) :: osobnik
767 character(80) :: lista
775 write(a,"(A6,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)") "WLONG=",os&
776 &obnik(1)," WSCP=",osobnik(2)," WELEC=",osobnik(3)," WBOND=",osobni&
777 &k(4)," WANG=",osobnik(5)," &"
779 write(a,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)") "WSCLOC=",o&
780 &sobnik(6)," WTOR=",osobnik(7)," WTORD=",osobnik(8)," WCORRH=",osob&
781 &nik(9),"WCORR5=",osobnik(11)," &"
783 write(a,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)") "WCORR6=",o&
784 &sobnik(12)," WEL_LOC=",osobnik(13)," WTURN3=",osobnik(14)," WTURN4&
785 &=",osobnik(15),"WTURN6=",osobnik(16)," &"
787 write(a,"(A7,F7.5,A6,F7.5,A8,F7.5,A)") "WVDWPP=",osobnik(17)," WHP&
788 &B=",osobnik(18)," WSCCOR=",osobnik(19)," &
791 write(a,"(A7,F7.5,A8,F7.5)") "CUTOFF=",coff," WCORR4=",osobnik(10)
796 c ======================================================================
797 c FunkcjaCeluR function
798 c ======================================================================
799 c Real scoring function
800 c ----------------------------------------------------------------------
802 real*8 function FunkcjaCeluR(zscore)
804 if (zscore.eq.0) then
807 FunkcjaCeluR = 1/zscore
812 c ======================================================================
813 c ReadInput subroutine
814 c ======================================================================
816 subroutine ReadInput(status)
823 character*100 :: wiersz,tmp
824 inquire(FILE=inpfn,EXIST=ex)
828 call write2log('==[ Reading main input ]======================')
832 read(inp, '(A)', iostat=stat) wiersz
835 if ((wiersz(1:4).eq.'PDB=').or.(wiersz(1:4).eq.'pdb=')) then
837 tmp = wiersz(5:len_trim(wiersz))
839 if (tmp(i:i).eq.' ') then
843 if (npdb.gt.maxnpdb) then
844 call write2log("Number of input PDB exceeds maxnpdb!")
849 if (index(trim(tmp),' ').gt.0) then
850 pdbfiles(i)=tmp(1:index(trim(tmp),' '))
852 pdbfiles(i)=tmp(1:len_trim(tmp))
854 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
856 endif ! Koniec czytania "PDB="
859 if ((wiersz(1:4).eq.'ALG=').or.(wiersz(1:4).eq.'alg=')) then
860 alg=wiersz(5:len_trim(wiersz))
863 call write2log ("Using 'simple' algorithm")
865 call write2log ("Using 'better' algorithm")
867 call write2log ("Using 'betterv2' algorithm")
869 call write2log ("Using 'CSA' algorithm")
871 call write2log ("Using 'Clustering CSA' algorithm")
874 call write2log ("Unknown algorithm. Using 'csa' as default")
879 select case (wiersz(1:12))
880 case('GENERATIONS=','generations=')
881 tmp = wiersz(13:len_trim(wiersz))
882 read(tmp,'(I)') maxgen
886 c select case(wiersz(1:9))
887 c case('CICUTOFF=','cicutoff=')
888 c tmp = wiersz(10:len_trim(wiersz))
889 c read(tmp(1:len_trim(tmp)),'(F7.5)') cicutoff
890 c call write2log("Initial CSA cutoff is set to "//tmp)
894 select case(wiersz(1:6))
895 case('MINCO=','minco=')
896 tmp = wiersz(7:len_trim(wiersz))
897 read(tmp(1:len_trim(tmp)),'(F7.5)') minco
898 call write2log("Minimal CSA cutoff factor is set to "//tmp)
902 select case(wiersz(1:6))
903 case('MAXCO=','maxco=')
904 tmp = wiersz(7:len_trim(wiersz))
905 read(tmp(1:len_trim(tmp)),'(F7.5)') maxco
906 call write2log("Maximal CSA cutoff factor is set to "//tmp)
911 select case(wiersz(1:11))
912 case('POPULATION=','population=')
913 tmp = wiersz(12:len_trim(wiersz))
914 read(tmp,'(I)') banksize
915 call write2log("Bank size is set to "//tmp)
919 select case(wiersz(1:13))
920 case('WHAMTEMPLATE=','whamtemplate=')
921 tmp = wiersz(14:len_trim(wiersz))
924 if (tmp(i:i).eq.' ') then
928 if (ntwham.gt.maxnpdb) then
929 call write2log("Number of WHAM templates exceeds maxnpdb!")
934 if (index(trim(tmp),' ').gt.0) then
935 whamtemplate(i)=tmp(1:index(trim(tmp),' '))
937 whamtemplate(i)=tmp(1:len_trim(tmp))
939 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
944 select case(wiersz(1:14))
945 case('MREMDTEMPLATE=','mremdtemplate=')
946 tmp = wiersz(15:len_trim(wiersz))
949 if (tmp(i:i).eq.' ') then
953 if (ntmremd.gt.maxnpdb) then
954 call write2log("Number of MREMD templates exceeds maxnpdb!")
959 if (index(trim(tmp),' ').gt.0) then
960 mremdtemplate(i)=tmp(1:index(trim(tmp),' '))
962 mremdtemplate(i)=tmp(1:len_trim(tmp))
964 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
969 select case(wiersz(1:7))
970 case('MAXMIN=','maxmin=')
971 tmp = wiersz(8:len_trim(wiersz))
972 read(tmp,'(I)') maxminstep
973 if (maxminstep.gt.0) then
974 call write2log("ZSCORE weights minimalization set to "//tmp)
977 call write2log("ZSCORE weights minimalization disabled")
983 select case(wiersz(1:8))
984 case('SCRIPTS=','scripts=')
986 tmp = wiersz(9:len_trim(wiersz))
988 if (tmp(i:i).eq.' ') then
992 if (nscripts.gt.maxscripts) then
993 call write2log("Number of scripts exceeds maxscripts!")
998 if (index(trim(tmp),' ').gt.0) then
999 scripts(i)=tmp(1:index(trim(tmp),' '))
1001 scripts(i)=tmp(1:len_trim(tmp))
1003 call write2log("Script will be used: "//scripts(i))
1004 tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
1006 end select ! Koniec czytania "scripts="
1011 c write additional info to log
1012 write(tmp,'(F8.5)') MUTRATIO
1013 call write2log("Mutation ratio is set to "//tmp)
1014 write(tmp,'(I2)') BANK_MULTIPLIER
1015 call write2log("Bank multiplier is set to "//tmp)
1016 write(tmp,'(I5)') BANK_MULTIPLIER*banksize
1018 call write2log("Trail population size will be "//trim(tmp)//" in&
1021 call write2log("Input file unresga.inp does not exist!!")
1022 write(*,*) "Input file unresga.inp does not exist!!"
1026 c Checking set parameters after reading input
1030 call write2log("Can't find 'PDB=' entry in input file")
1034 inquire(FILE=trim(pdbfiles(i)),EXIST=ex)
1036 call write2log("Can't find file "//trim(pdbfiles(i)))
1043 if (maxgen.eq.0) then
1044 call write2log("Couldn't find GENERATIONS= entry in input file")
1049 if (banksize.eq.0) then
1050 call write2log("Can't find POPULATION= entry in input file")
1052 elseif(mod(banksize,2).eq.1) then
1053 call write2log("POPULATION has to be a even number.")
1058 if (ntwham.ne.npdb) then
1059 call write2log("Number of WHAM templates and PDB files is not equ&
1065 if (ntmremd.ne.npdb) then
1066 call write2log("Number of MREMD templates and PDB files is not eq&
1072 if (nscripts.gt.0) then
1074 inquire(FILE=trim(scripts(i)),EXIST=ex)
1076 call write2log("Can't find file "//trim(scripts(i)))
1083 end subroutine ReadInput
1086 c ======================================================================
1087 c CreateInputs subroutine
1088 c ======================================================================
1089 c Creates input files for MREMD, WHAM and ZSCORE.
1090 c ----------------------------------------------------------------------
1092 subroutine CreateInputs(rozmiar,pop)
1095 include 'common.inc'
1096 character(80),dimension(5) :: wagi
1097 integer*4 :: rozmiar,i,j,k
1099 dimension pop(rozmiar,21)
1100 !character*50 :: pdbfiles(maxnpdb)
1101 character*50 :: command,tmptext, prefix
1102 character*256 :: line,plik
1103 integer :: gdzie = 0
1104 character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1105 &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&
1107 !write(tmptext,'(I10)') banksize
1108 !call write2log(trim(tmptext))
1111 prefix=pdbfiles(j)(1:len_trim(pdbfiles(j))-4)
1112 c Create MREMD input
1113 open(imremd, file = mremdtemplate(j))
1114 plik=trim(prefix)//"_"//omremdfn
1115 call write2log(trim("Creating MREMD input: "//plik))
1116 open(omremd, file = plik)
1117 3000 read(imremd, '(A)', end = 3010) line
1118 gdzie = index (line, '{PREFIX}')
1119 if (gdzie .ne. 0) then
1120 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1121 write(omremd,"(A)") trim(line)
1122 elseif (index(line,'{POPSIZE}') .ne. 0) then
1123 gdzie = index (line, '{POPSIZE}')
1124 write(tmptext,'(I4)') rozmiar
1125 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(
1127 write(omremd,"(A)") trim(line)
1128 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1129 gdzie = index (line, '{WEIGHTS}')
1131 call Weights2Str(pop(i,:),CUTOFF,wagi)
1132 write(omremd,"(A)") wagi(1)
1133 write(omremd,"(A)") wagi(2)
1134 write(omremd,"(A)") wagi(3)
1135 write(omremd,"(A)") wagi(4)
1136 write(omremd,"(A)") wagi(5)
1139 write(omremd,"(A)") trim(line)
1148 open(iwham, file = whamtemplate(j))
1149 plik=trim(prefix)//"_"//owhamfn
1150 call write2log(trim("Creating WHAM input: "//plik))
1151 open(owham, file = plik)
1152 3015 read(iwham, '(A)', end = 3020) line
1153 gdzie = index (line, '{PREFIX}')
1154 if (gdzie .ne. 0) then
1155 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
1156 write(owham,"(A)") trim(line)
1157 elseif (index(line,'{POPSIZE}') .ne. 0) then
1158 gdzie = index (line, '{POPSIZE}')
1159 write(tmptext,'(I4)') rozmiar
1160 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1162 write(owham,"(A)") trim(line)
1163 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1164 gdzie = index (line, '{WEIGHTS}')
1166 call Weights2Str(pop(i,:),CUTOFF,wagi)
1167 write(owham,"(A)") wagi(1)
1168 write(owham,"(A)") wagi(2)
1169 write(owham,"(A)") wagi(3)
1170 write(owham,"(A)") wagi(4)
1171 write(owham,"(A)") wagi(5)
1172 write(owham,"(A)") ''
1175 write(owham,"(A)") trim(line)
1183 c Create ZSCORE input
1184 open(izscore, file = izscorefn)
1185 call write2log(trim("Creating ZSCORE input: "//ozscorefn))
1186 open(ozscore, file = ozscorefn)
1187 3025 read(izscore, '(A)', end = 3030) line
1188 gdzie = index (line, '{PREFIX')
1189 if (gdzie .ne. 0) then
1190 read(line(gdzie+7:gdzie+9),'(I2.2)') k
1191 prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1192 line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+10:)&
1195 write(ozscore,"(A)") trim(line)
1196 elseif (index(line,'{PARAM') .ne. 0) then
1197 gdzie=index(line,'{PARAM')
1198 read(line(gdzie+6:gdzie+8),'(I2.2)') k
1199 prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
1201 write(tmptext,'(I3.3)') i
1202 line=trim(line(1:gdzie-1))//trim(prefix)//"_wham_par"//tmptex&
1204 write(ozscore,"(A)") trim(line)
1206 elseif (index(line,'{POPSIZE}') .ne. 0) then
1207 gdzie = index (line, '{POPSIZE}')
1208 write(tmptext,'(I4)') rozmiar
1209 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1211 write(ozscore,"(A)") trim(line)
1212 elseif (index(line,'{WEIGHTS}') .ne. 0) then
1213 gdzie = index (line, '{WEIGHTS}')
1215 call Weights2Str(pop(i,:),CUTOFF,wagi)
1216 write(ozscore,"(A)") wagi(1)
1217 write(ozscore,"(A)") wagi(2)
1218 write(ozscore,"(A)") wagi(3)
1219 write(ozscore,"(A)") wagi(4)
1220 write(ozscore,"(A)") trim(wagi(5))
1221 write(ozscore,"(A)") ''
1223 elseif (index(line,'{MAXMIN}') .ne. 0) then
1224 gdzie = index (line, '{MAXMIN}')
1225 write(tmptext,'(I4)') maxminstep
1226 line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
1228 write(ozscore,"(A)") trim(line)
1230 write(ozscore,"(A)") trim(line)
1239 call write2log('Copying scripts is disabled.')
1241 c ! Open output file
1242 c plik=trim(prefix)//"/"//trim(scripts(i))
1243 c open(oscript, file = plik)
1245 c plik=trim(scripts(i))
1246 c open(iscript, file = plik)
1247 c ! rewrite from input to output
1248 c3030 read(iscript, '(A)', end = 3035) line
1249 c write(oscript,'(A)') trim(line)
1257 c ==========================================
1259 c ==========================================
1262 c write(ow,'(I4)') generation
1263 c write(ow,'(200A)') "# WLONG WSCP WELEC WBOND WANG WSC&
1264 c &LOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC WTURN&
1265 c &3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR"
1267 c write(ow,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),pop&
1268 c &(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),pop&
1269 c &(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,19)
1273 end subroutine CreateInputs
1276 c ======================================================================
1277 c ZnajdzPodobnego function
1278 c ======================================================================
1279 c Finds the most similar individual to osobnik in population(populacja)
1280 c and returns his position in the population
1281 c ----------------------------------------------------------------------
1283 integer*4 function ZnajdzPodobnego(rozmiar,populacja,osobnik,csacu&
1285 integer*4 :: rozmiar,pozycja
1286 real*8 :: delta,mindelta,csacutoff
1288 real*8, dimension(rozmiar,21) :: populacja
1289 real*8, dimension(21) :: osobnik
1291 mindelta = csacutoff
1292 c mindelta = 10000.0
1298 delta=delta+(populacja(i,j)-osobnik(j))**2
1300 c write(tmp,'(2F7.5 )') delta,csacutoff
1301 c call write2log("Delta is "//tmp)
1303 if (delta.lt.mindelta) then
1306 write(tmp,'(F7.5,X,I)') delta,pozycja
1307 call write2log("Delta is "//tmp)
1311 ZnajdzPodobnego = pozycja
1314 c ======================================================================
1315 c WybierzOsobnika function
1316 c ======================================================================
1317 c Selects individual from population using roulette method
1318 c ----------------------------------------------------------------------
1320 integer*4 function WybierzOsobnika(rozmiar,populacja)
1321 integer*4 :: osobnik, rozmiar
1322 real*8 :: los, partsum
1324 real*8, dimension(rozmiar,21) :: populacja
1330 osobnik = osobnik + 1
1331 partsum = partsum + populacja(osobnik,21)
1332 if ((osobnik.eq.rozmiar).or.(partsum.ge.los)) exit
1335 WybierzOsobnika = osobnik
1339 c ======================================================================
1340 c Najlepszy function
1341 c ======================================================================
1342 c Selects the individual with the best score from population
1343 c ----------------------------------------------------------------------
1345 integer*4 function Najlepszy(rozmiar, maks, populacja)
1346 integer*4 :: rozmiar
1348 real*8, dimension(rozmiar,21) :: populacja
1351 if (populacja(i,20).eq.maks) then
1358 c ======================================================================
1359 c WriteState subroutine
1360 c ======================================================================
1361 c Writes current state of calculations in case jobs have to be
1363 c ----------------------------------------------------------------------
1365 subroutine WriteState()
1368 include 'common.inc'
1370 open(ostate,file=ostatefn)
1371 write(ostate,'(I4)') generation
1372 write(ostate,'(L2)') do_optima
1373 write(ostate,'(L2)') do_ga
1374 write(ostate,'(F7.5)') avrd
1376 end subroutine WriteState
1378 c ======================================================================
1379 c ReadPopulation subroutine
1380 c ======================================================================
1381 c Reads population from population.summary file
1382 c ----------------------------------------------------------------------
1384 subroutine ReadPopulation(nind,pop)
1387 character*8 :: dummy
1388 dimension pop(nind,21)
1391 character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
1392 &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&
1394 read(opopsum,'(A)') dummy
1396 read(opopsum,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),p&
1397 &op(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),p&
1398 &op(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,1&
1401 end subroutine ReadPopulation
1403 c ======================================================================
1404 c ReadOptimaW subroutine
1405 c ======================================================================
1406 c Reads weights from zscore generated weights_opt.zscorepar_XXX files
1408 c nind - number of individuals/files to read
1409 c pop - array where weights are stored
1410 c ----------------------------------------------------------------------
1412 subroutine ReadOptimaW(nind,pop)
1413 integer*4 :: i,nind,stat
1416 character*40 :: filename
1417 character*200 :: tmptxt
1418 dimension pop(nind,21)
1423 c Cosmetic stuff for pretty log
1424 call write2log("==[ Reading weights from zscore files ]=======")
1425 call write2log("# WLONG WSCP WELEC WBOND WANG WS&
1426 &CLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_L&
1427 &OC WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR")
1429 c Actual weight read and error handling
1432 write(tmptxt,'(I3.3)') i
1433 filename = "weights_opt.zscorepar_"//trim(tmptxt)
1434 open(izsoptw, file = filename)
1435 read(izsoptw,"(A4,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)",iosta&
1436 &t=stat) tmptxt,pop(i,1),tmptxt,pop(i,2),tmptxt,pop(i,3),tmptxt,pop&
1437 &(i,4),tmptxt,pop(i,5)
1440 call write2log("-- ERROR while reading "//trim(filename)//"--")
1446 read(izsoptw,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1447 &t=stat) tmptxt,pop(i,6),tmptxt,pop(i,7),tmptxt,pop(i,8),tmptxt,pop&
1448 &(i,9),tmptxt,pop(i,11)
1451 call write2log("-- ERROR while reading "//trim(filename)//"--")
1457 read(izsoptw,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
1458 &t=stat) tmptxt,pop(i,12),tmptxt,pop(i,13),tmptxt,pop(i,14),tmptxt,&
1459 &pop(i,15),tmptxt,pop(i,16)
1462 call write2log("-- ERROR while reading "//trim(filename)//"--")
1468 read(izsoptw,"(A7,F7.5,A8,F7.5,A6,F7.5,A)",iostat=stat) tmptxt, &
1469 &pop(i,19),tmptxt,pop(i,17),tmptxt,pop(i,18),tmptxt
1472 call write2log("-- ERROR while reading "//trim(filename)//"--")
1478 read(izsoptw,"(A22,F7.5)",iostat=stat) tmptxt, pop(i,10)
1481 call write2log("-- ERROR while reading "//trim(filename)//"--")
1488 write(tmptxt,'(19F9.5)') (pop(i,j),j=1,19)
1489 call write2log(trim(tmptxt))
1494 end subroutine ReadOptimaW
1497 c ======================================================================
1498 c ReadZEnergy subroutine
1499 c ======================================================================
1500 c Reads Energy from ZSCORE output files zscore_out_parXXX_GB000
1501 c ----------------------------------------------------------------------
1503 subroutine ReadZEnergy(nind,pop)
1507 character*40 :: filename,tmptxt,tmptxt2
1508 character*100 :: wiersz
1509 dimension pop(nind,21)
1512 call write2log("==[ Reading Final Function Value from zscore ]===&
1515 write(tmptxt,'(I3.3)') i
1517 filename=trim(izenergyfn)//trim(tmptxt)//"_GB000"
1518 open(izenergy, file = filename)
1520 read(izenergy, '(A)', iostat=stat) wiersz
1522 if (wiersz(1:22).eq.' Final function value:') then
1523 tmptxt=trim(adjustl(wiersz(23:48)))
1524 read(tmptxt,'(F16.5)') pop(i,20)
1525 write(tmptxt2,'(E15.7)') pop(i,20)
1526 call write2log("Reading score from "//filename//" : "//tmptxt/&
1530 if (pop(i,20).eq.0) then
1531 c Setting FFV to a realy big valiue or otherwise CalcFittness will have problems
1532 pop(i,20)=huge(0.0d0)
1533 call write2log("ERROR while reading FFV from "//filename)
1538 end subroutine ReadZEnergy
1540 c =======================================================================
1541 c WritePopSum subroutine
1542 c =======================================================================
1543 c Writes whole individuals (weights, score, ID, generation) to file
1544 c population.summary This file can be used to plot evolution of the
1545 c weights and scores over time.
1546 c -----------------------------------------------------------------------
1547 subroutine WritePopSum()
1548 logical :: jestplik = .false.
1549 character*5 :: tmptext
1552 c dimension bank(nind,21)
1554 include 'common.inc'
1557 c write(*,*) banksize
1559 write(tmptext,'(I3.3)') generation
1560 inquire(file=opopsumfn, exist=jestplik)
1561 if (.not. jestplik) then
1562 open(opopsum, file = opopsumfn)
1563 write(opopsum, "(A189)") "# WLONG WSCP WELEC WBOND WANG &
1564 & WSCLOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 WEL_LOC&
1565 & WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCOR SCO&
1570 open(opopsum, file = opopsumfn, access = "append")
1571 c call ReadPopulation(banksize,populacja)
1572 write(opopsum, "(A189)") "# ==== Generation "//trim(tmptext)//" ==&
1573 &==================================================================&
1574 &==================================================================&
1575 &================================="
1576 c call write2log(banksize)
1579 c write(opopsum,"(F9.5,2F8.5)") populacja(i,1),populacja(i,2),popu&
1581 write(opopsum,"(F9.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,&
1582 &F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F20.5,I8,I8)") bank(i&
1583 &,1),bank(i,2),bank(i,3),bank(i,4),bank(i,5),bank(i,6),bank(i,7),ba&
1584 &nk(i,8),bank(i,9),bank(i,10),bank(i,11),bank(i,12),bank(i,13),bank&
1585 &(i,14),bank(i,15),bank(i,16),bank(i,17),bank(i,18),bank(i,19),bank&
1591 end subroutine WritePopSum
1593 c ======================================================================
1594 c WriteBank subroutine
1595 c ======================================================================
1596 c Writes CSA BANK to file
1597 c ----------------------------------------------------------------------
1598 subroutine WriteBank(b)
1600 include 'common.inc'
1601 real*8,dimension(banksize,21) :: b
1602 character*250 :: tmptext
1603 character*250 :: header
1606 header="# WLONG WSCP WELEC WBOND WANG WSC&
1607 &LOC WTOR WTORD WCORRH WCORR4 WCORR5 WCORR6 &
1608 &WEL_LOC WTURN3 WTURN4 WTURN6 WVDWPP WHPB WSCCO&
1611 call write2log("Writing Bank to file "//iobankfn)
1612 open(iobank, file = iobankfn)
1613 write(iobank, "(A)") trim(adjustl(header))
1616 write(iobank,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1617 write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1618 call write2log(trim(tmptext))
1623 c ======================================================================
1624 c ReadBank subroutine
1625 c ======================================================================
1626 c Read CSA BANK from file
1627 c ----------------------------------------------------------------------
1628 subroutine ReadBank(b)
1630 include 'common.inc'
1631 real*8,dimension(banksize,21) :: b
1632 character*250 :: tmptext
1634 call write2log("Reading Bank from file "//iobankfn)
1635 open(iobank, file = iobankfn)
1637 read(iobank,"(A)") tmptext
1639 read(iobank,"(19F10.5,E15.7,F10.5)") (b(i,j),j=1,21)
1640 write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
1641 call write2log(trim(tmptext))
1646 c ======================================================================
1647 c CalcFitness subroutine
1648 c ======================================================================
1649 c Calculates fitness of every individual in the bank
1650 c ----------------------------------------------------------------------
1651 subroutine CalcFitness(b,fitn)
1652 include 'common.inc'
1653 real*8,dimension(banksize,21) :: b
1654 real*8 :: fitn, sumfitn
1660 fitn=fitn+(1/b(i,20))
1661 sumfitn=sumfitn+b(i,20)
1664 b(i,21)=1/(b(i,20)*fitn)
1665 c b(i,21)=b(i,20)/fitn
1670 c ======================================================================
1671 c CalcAvgDist subroutine
1672 c ======================================================================
1673 c Calculates average distance between individuals in the bank
1674 c ----------------------------------------------------------------------
1675 subroutine CalcAvgDist(b,avgd)
1676 include 'common.inc'
1677 real*8,dimension(banksize,21) :: b
1682 nd = (banksize-1)*banksize/2 ! number of distances to calculate
1687 d=d+(b(i,w)-b(j,w))**2
1694 c -----------------------------------------------------------------------