1 subroutine read_general_data(*)
4 include "DIMENSIONS.ZSCOPT"
8 integer ierror,kolor,klucz
10 include "COMMON.WEIGHTS"
11 include "COMMON.NAMES"
12 include "COMMON.VMCPAR"
13 include "COMMON.TORSION"
14 include "COMMON.INTERACT"
15 include "COMMON.ENERGIES"
16 include "COMMON.MINPAR"
17 include "COMMON.IOUNITS"
18 include "COMMON.TIME1"
19 include "COMMON.PROTFILES"
20 include "COMMON.CHAIN"
21 include "COMMON.CLASSES"
22 include "COMMON.THERMAL"
23 include "COMMON.FFIELD"
24 include "COMMON.OPTIM"
25 include "COMMON.CONTROL"
26 include "COMMON.SCCOR"
27 include "COMMON.SPLITELE"
28 character*800 controlcard
29 integer i,j,k,l,ii,n_ene_found,ist1,iet1,ist2,iet2,ls,le
30 integer ind,itype1,itype2,itypf,itypsc,itypp,itypt1,itypt2
31 integer ilen,lenpot,lenpre
33 character*3 liczba,liczba1
40 double precision xchuj,weitemp,weitemp_low,weitemp_up
45 c write (iout,*) "Enter read_general_data"
53 C Read first record: seed and number of energy components
54 call card_concat(controlcard,.true.)
55 c write (iout,*) "card1",controlcard
57 call readi(controlcard,"ISEED",iseed,0)
58 if (iseed.eq.0) stop "Seed is zero!"
59 c print *,me," iseed",iseed
60 call readi(controlcard,"NPARMSET",nparmset,1)
61 c print *,me," nparmset",nparmset
63 c Split processor pool if multiple parameter sets are treated
64 if (nparmset.eq.1) then
66 c print *,me," opening ",outname," opened"
67 if (me.eq.Master)open(iout,file=outname,status='unknown')
69 c print *,me," outname ",outname," opened"
71 if (nprocs.lt.nparmset) then
72 if (me.eq.Master) write (*,*)
73 & "*** Cannot split parameter sets for fewer processors than sets",
75 call MPI_Finalize(ierror)
78 c write (iout,*) "nparmset",nparmset
79 nprocs = nprocs/nparmset
81 klucz = mod(me,nprocs)
82 c write (*,*) "My old rank",me," kolor",kolor," klucz",klucz
83 call MPI_Comm_split(ALL_COMM,kolor,klucz,WHAM_COMM,ierror)
84 call MPI_Comm_size(WHAM_COMM,nprocs,ierror)
85 call MPI_Comm_rank(WHAM_COMM,me,ierror)
86 c write (*,*) "My new rank",me," comm size",nprocs
87 c write (*,*) "ALL_COMM",ALL_COMM
88 c & " WHAM_COMM",WHAM_COMM
90 c write (*,*) "My parameter set is",myparm
91 write(liczba,'(bz,i3.3)') me
92 write(liczba1,'(bz,i3.3)') myparm
93 outname=prefix(:lenpre)//'.out_par'//liczba1//'_'//
94 & pot(:lenpot)//liczba
95 if (me.eq.Master)open(iout,file=outname,status='unknown')
98 c print *,me," checkpoint 1"
99 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
100 call readi(controlcard,'TORMODE',tor_mode,0)
101 c print *,me," tor_mode",tor_mode
102 call readi(controlcard,"N_ENE",n_ene,max_ene)
103 c print *,"iseed",iseed," n_ene",n_ene
104 call readi(controlcard,"NPARMSET",nparmset,1)
105 geom_and_wham_weights =
106 & index(controlcard,"GEOM_AND_WHAM_WEIGHTS").gt.0
107 c write (iout,*) "GEOM_AND_WHAM_WEIGHTS",geom_and_wham_weights
108 if (iseed.gt.0) iseed=-iseed
110 out_newe=index(controlcard,"OUT_NEWE").gt.0
111 unres_pdb = index(controlcard,"UNRES_PDB").gt.0
112 c write (iout,*) "UNRES_PDB ",unres_pdb
113 c Energy calculation settings section
114 call readi(controlcard,'TORMODE',tor_mode,0)
115 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
116 call reada(controlcard,'BOXX',boxxsize,100.0d0)
117 call reada(controlcard,'BOXY',boxysize,100.0d0)
118 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
119 c print *,me," checkpoint 2"
120 c Cutoff range for interactions
121 call reada(controlcard,"R_CUT",r_cut,15.0d0)
122 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
123 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
124 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
125 if (lipthick.gt.0.0d0) then
126 bordliptop=(boxzsize+lipthick)/2.0
127 bordlipbot=bordliptop-lipthick
129 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)
131 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
132 buflipbot=bordlipbot+lipbufthick
133 bufliptop=bordliptop-lipbufthick
134 if ((lipbufthick*2.0d0).gt.lipthick.and.me.eq.Master)
135 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
137 c write(iout,*) "bordliptop=",bordliptop
138 c write(iout,*) "bordlipbot=",bordlipbot
139 c write(iout,*) "bufliptop=",bufliptop
140 c write(iout,*) "buflipbot=",buflipbot
141 call readi(controlcard,'SYM',symetr,1)
142 c print *,me," checkpoint 3"
143 c write (iout,*) "n_ene",n_ene
145 c Energy-term-weights section
147 C Read third record: weights
151 c print *,me," checkpoint 4"
152 C Read fourth record: masks
153 call card_concat(controlcard,.true.)
154 c write (iout,*) "card2",controlcard
156 key = "MASK_"//wname(i)(2:ilen(wname(i)))
157 call readi(controlcard,key(:ilen(key)),imask(i),0)
159 c print *,me," checkpoint 5"
160 C Read fifth record: lower limits of weights
161 call card_concat(controlcard,.true.)
162 c write (iout,*) "card3",controlcard
164 key = "WLOW_"//wname(i)(2:ilen(wname(i)))
165 call reada(controlcard,key(:ilen(key)),ww_low(i),ww(i))
167 C Read sixth record: upper limits of weights
168 call card_concat(controlcard,.true.)
169 c write (iout,*) "card4",controlcard
171 key = "WUP_"//wname(i)(2:ilen(wname(i)))
172 call reada(controlcard,key(:ilen(key)),ww_up(i),ww(i))
174 C Read seventh record: VMC parameters
175 call card_concat(controlcard,.true.)
176 c write (iout,*) "card5",controlcard
177 call readi(controlcard,"MODE",mode,3)
178 call readi(controlcard,"MAXSTEP_SCAN",maxstep_scan,50)
179 call reada(controlcard,"RSTEP_SCAN",step_scan,1.0d-1)
180 call readi(controlcard,"READ_STAT",read_stat,3)
181 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
182 c init_ene = index(controlcard,"INIT_ENE").gt.0 .and. read_stat.gt.1
184 call readi(controlcard,"NMCM",nmcm,0)
185 call readi(controlcard,"MAXSCAN",maxscan,0)
186 call readi(controlcard,"MAXMIN",maxmin,100)
187 call readi(controlcard,"MAXFUN",maxfun,100)
188 call reada(controlcard,"TOLF",tolf,1.0d-6)
189 call reada(controlcard,"RTOLF",rtolf,1.0d-6)
191 if (index(controlcard,"OUT_MINIM").gt.0) out_minim=iout
193 if (index(controlcard,"PRINT_INI").gt.0) print_ini=1
195 if (index(controlcard,"PRINT_FIN").gt.0) print_fin=1
197 if (index(controlcard,"PRINT_STAT").gt.0) print_stat=1
198 call reada(controlcard,"RSTIME",rstime,9.0d2)
199 call reads(controlcard,"MINIMIZER",minimizer,"SUMSL")
200 call readi(controlcard,"OPT_MODE",opt_mode,0)
201 mod_other_params=index(controlcard,"OPTIMIZE_OTHER").gt.0
202 if (read_stat.eq.0 .and. mod_other_params) then
204 & write (iout,*) "Error: only optimization of energy-term",
205 & " weights can be performed with READ_STAT=",read_stat
209 if (index(controlcard,"RESTART").gt.0) then
214 c print *,me," checkpoint 6"
217 c-----------------------------------------------------------------------
218 subroutine read_optim_parm(*)
221 include "DIMENSIONS.ZSCOPT"
225 integer ierror,kolor,klucz
227 include "COMMON.WEIGHTS"
228 include "COMMON.NAMES"
229 include "COMMON.VMCPAR"
230 include "COMMON.TORSION"
231 include "COMMON.INTERACT"
232 include "COMMON.ENERGIES"
233 include "COMMON.MINPAR"
234 include "COMMON.IOUNITS"
235 include "COMMON.TIME1"
236 include "COMMON.PROTFILES"
237 include "COMMON.CHAIN"
238 include "COMMON.CLASSES"
239 include "COMMON.THERMAL"
240 include "COMMON.FFIELD"
241 include "COMMON.OPTIM"
242 include "COMMON.CONTROL"
243 include "COMMON.SCCOR"
244 character*800 controlcard
246 integer i,j,k,l,ii,n_ene_found,ist1,iet1,ist2,iet2,ls,le
247 integer ind,ind1,ind2,itype1,itype2,itypf,itypsc,itypp,
248 & itypt1,itypt2,masktemp,iblock,iaux
249 integer ilen,lenpot,lenpre
252 character*3 liczba,liczba1
259 double precision xchuj,weitemp,weitemp_low,weitemp_up
268 & write (iout,*) "MOD_OTHER_PARAMS ",mod_other_params
272 mask_tor(0,itypt1,itypt2,iblock)=0
273 weitor(0,itypt1,itypt2,iblock)=1.0d0
274 weitor_low(0,itypt1,itypt2,iblock)=1.0d0
275 weitor_up(0,itypt1,itypt2,iblock)=1.0d0
281 mask_tor(l,itypt1,itypt2,iblock)=0
282 weitor(l,itypt1,itypt2,iblock)=1.0
283 weitor_low(l,itypt1,itypt2,iblock)=1.0
284 weitor_up(l,itypt1,itypt2,iblock)=1.0
290 write (iout,*) "ntyp",ntyp
293 write (iout,*) "itypt1",itypt1," itypt2",itypt2,
294 & "weitor",weitor(0,itypt1,itypt2,1),eitor(0,itypt1,itypt2,2)
298 if (mod_other_params) then
300 c &"Internal parameters of UNRES energy components will be optimized"
301 call card_concat(controlcard,.true.)
302 c write (iout,*) "mod_side ",controlcard
304 mod_side = (index(controlcard,"MOD_SIDE").gt.0)
308 call card_concat(controlcard,.true.)
309 do while ( index(controlcard,"END").eq.0 )
310 call split_string(controlcard,item(1),4,nitem)
311 if (nitem.eq.1 .or. item(2)(:1).eq."*") then
312 nsingle_sc=nsingle_sc+1
313 ityp_ssc(nsingle_sc)=rescode(1,item(1),0)
315 epss_low(nsingle_sc)=0.02d0
317 read (item(3),*) epss_low(nsingle_sc)
320 epss_up(nsingle_sc)=0.0d0
322 read (item(4),*) epss_up(nsingle_sc)
326 ityp_psc(1,npair_sc)=rescode(1,item(1),0)
327 ityp_psc(2,npair_sc)=rescode(1,item(2),0)
329 epsp_low(npair_sc)=0.02d0
331 read (item(3),*) epsp_low(npair_sc)
334 epsp_up(npair_sc)=0.0d0
336 read (item(4),*) epsp_up(npair_sc)
339 call card_concat(controlcard,.true.)
341 if (nsingle_sc+npair_sc.eq.0) mod_side=.false.
342 call card_concat(controlcard,.true.)
345 & index(controlcard,"MOD_SIDE_OTHER").gt.0
346 c write (iout,*) "mod_side_other ",controlcard,mod_side_other
347 if (mod_side_other) then
348 mod_side_other=.false.
354 call card_concat(controlcard,.true.)
355 c write (iout,*) "mod_side_oter ",controlcard
356 do while ( index(controlcard,"END").eq.0 )
357 call reads(controlcard,"RESKIND",reskind," ")
358 itypsc=rescode(1,reskind,0)
359 if (itypsc.lt.1 .or. itypsc.gt.20) then
360 if (me.eq.Master) then
361 write (iout,*) "Error in SC optimization data;",
362 & " SC type must not be dummy, type is" ,restyp," ",itypsc
363 write (*,*) "Error in SC optimization data;",
364 & " SC type must not be dummy, type is" ,restyp," ",itypsc
368 call readi(controlcard,"MASK_SIGMA0",mask_sigma(1,itypsc),0)
369 call readi(controlcard,"MASK_SIGMAII",mask_sigma(2,itypsc),0)
370 call readi(controlcard,"MASK_CHIP",mask_sigma(3,itypsc),0)
371 call readi(controlcard,"MASK_ALP",mask_sigma(4,itypsc),0)
372 call reada(controlcard,"SIGMA0_LOW",sigma_low(1,itypsc),0.d0)
373 call reada(controlcard,"SIGMA0_UP",sigma_up(1,itypsc),0.0d0)
374 call reada(controlcard,"SIGMAII_LOW",sigma_low(2,itypsc),
376 call reada(controlcard,"SIGMAII_UP",sigma_up(2,itypsc),0.0d0)
377 call reada(controlcard,"CHIP_LOW",sigma_low(3,itypsc),0.d0)
378 call reada(controlcard,"CHIP_UP",sigma_up(3,itypsc),0.0d0)
379 call reada(controlcard,"ALP_LOW",sigma_low(4,itypsc),0.d0)
380 call reada(controlcard,"ALP_UP",sigma_up(4,itypsc),0.0d0)
382 if (sigma_low(k,itypsc).eq.0.0d0 .and.
383 & sigma_up(k,itypsc).eq.0.0d0) mask_sigma(k,itypsc)=0
386 mod_side_other=mod_side_other.or.mask_sigma(k,itypsc).gt.0
389 & write (iout,'(a4,i3,4(i2,2f8.3))') reskind,itypsc,
390 & (mask_sigma(k,itypsc),
391 & sigma_low(k,itypsc),sigma_up(k,itypsc),k=1,4)
392 call card_concat(controlcard,.true.)
393 c write (iout,*) "mod_side_oter ",controlcard
395 if (me.eq.Master) then
396 write (iout,*) "Optimization flags of one-body SC parameters"
398 write (iout,'(a4,i3,4(i2,2f8.3))') restyp(i),i,
399 & (mask_sigma(k,i),sigma_low(k,i),sigma_up(k,i),k=1,4)
402 call card_concat(controlcard,.true.)
404 c write (iout,*) "mod_side_other ",mod_side_other
405 c write (iout,*) "mod_tor ",controlcard
407 mod_tor = index(controlcard,"MOD_TOR").gt.0
410 do i=-ntortyp,ntortyp
411 do j=-ntortyp,ntortyp
412 mask_tor(0,i,j,iblock)=0
413 weitor(0,i,j,iblock)=1.0d0
414 weitor_low(0,i,j,iblock)=0.0d0
415 weitor_up(0,i,j,iblock)=2.0d0
419 call card_concat(controlcard,.true.)
420 c write (iout,*) controlcard
421 do while ( index(controlcard,"END").eq.0 )
422 call split_string(controlcard,item(1),7,nitem)
425 & write (*,*) "Error in torsional optimization data;",
426 & " must specify both residues and type"
434 c write (iout,*) "item3 ",item(3)," item4 ",item(4),
437 if (nitem.gt.2) read(item(3),*) masktemp
438 if (nitem.gt.3) read(item(4),*) weitemp
439 if (nitem.gt.4) read(item(5),*) weitemp_low
440 if (nitem.gt.5) read(item(6),*) weitemp_up
441 if (nitem.gt.6) read(item(7),*) iblock
442 c write (iout,*) controlcard
443 c write (iout,*) item(1)," ",item(2),weitemp,
444 c & weitemp_low,weitemp_up
445 if (index(item(1),"*").gt.0) then
449 ist1=itortyp(rescode(1,item(1),0))
452 if (index(item(2),"*").gt.0) then
456 ist2=itortyp(rescode(1,item(2),0))
459 c write (iout,*) "ist1",ist1," iet1",iet1," ist2",ist2,
464 c write (iout,*) "itypt1",itypt1," itypt2",itypt2,
465 c & weitemp,weitemp_low,weitemp_up
466 mask_tor(0,itypt1,itypt2,iblock)=masktemp
467 weitor(0,itypt1,itypt2,iblock)=weitemp
468 weitor_low(0,itypt1,itypt2,iblock)=weitemp_low
469 weitor_up(0,itypt1,itypt2,iblock)=weitemp_up
470 c write (iout,*) "itypt1",itypt1," itypt2",itypt2,
471 c & mask_tor(0,itypt1,itypt2,iblock),
472 c & weitor(0,itypt1,itypt2,iblock),
473 c & weitor_low(0,itypt1,itypt2,iblock),
474 c & weitor_up(0,itypt1,itypt2,iblock)
477 call card_concat(controlcard,.true.)
478 c write (iout,*) controlcard
481 if (tor_mode.gt.1) then
482 if (me.eq.Master) write (iout,*) "TORMODE is",tor_mode,
483 & " torsional constants are computed from energy map expansion."
488 & write (iout,*) "Initialized torsional parameters:"
490 do itypt1=-ntortyp,ntortyp
491 do itypt2=-ntortyp,ntortyp
492 if (me.eq.Master)write(iout,'(3i5,3f10.5)') itypt1,itypt2,
493 & mask_tor(0,itypt1,itypt2,iblock),
494 & weitor(0,itypt1,itypt2,iblock),
495 & weitor_low(0,itypt1,itypt2,iblock),
496 & weitor_up(0,itypt1,itypt2,iblock)
501 if (tor_mode.eq.1) then
502 c Exchange indices because the residue order in new torsionals is reversed
504 do itypt1=-ntortyp,ntortyp
505 do itypt2=itypt1+1,ntortyp
506 iaux=mask_tor(0,itypt1,itypt2,iblock)
507 mask_tor(0,itypt1,itypt2,iblock)=
508 & mask_tor(0,itypt2,itypt1,iblock)
509 mask_tor(0,itypt2,itypt1,iblock)=iaux
510 aux=weitor(0,itypt1,itypt2,iblock)
511 weitor(0,itypt1,itypt2,iblock)=
512 & weitor(0,itypt2,itypt1,iblock)
513 weitor(0,itypt2,itypt1,iblock)=aux
514 aux=weitor_low(0,itypt1,itypt2,iblock)
515 weitor_low(0,itypt1,itypt2,iblock)=
516 & weitor_low(0,itypt2,itypt1,iblock)
517 weitor_low(0,itypt2,itypt1,iblock)=aux
518 aux=weitor_up(0,itypt1,itypt2,iblock)
519 weitor_up(0,itypt1,itypt2,iblock)=
520 & weitor_up(0,itypt2,itypt1,iblock)
521 weitor_up(0,itypt2,itypt1,iblock)=aux
526 call card_concat(controlcard,.true.)
528 c write (iout,*) "mod_sccor ",controlcard
530 mod_sccor = index(controlcard,"MOD_SCCOR").gt.0
532 call card_concat(controlcard,.true.)
535 do i=-nsccortyp,nsccortyp
536 do j=-nsccortyp,nsccortyp
537 mask_tor(l,i,j,iblock)=0
538 weitor(l,i,j,iblock)=1.0d0
539 weitor_low(l,i,j,iblock)=0.0d0
540 weitor_up(l,i,j,iblock)=2.0d0
545 do while ( index(controlcard,"END").eq.0 )
546 call split_string(controlcard,item(1),7,nitem)
549 & write (*,*) "Error in torsional optimization data;",
550 & " must specify both residues and type"
556 if (nitem.gt.3) read(item(4),*) masktemp
557 if (nitem.gt.4) read(item(5),*) weitemp
558 if (nitem.gt.5) read(item(6),*) weitemp_low
559 if (nitem.gt.6) read(item(7),*) weitemp_up
560 if (index(item(1),"*").gt.0) then
564 ist1=isccortyp(rescode(1,item(1),0))
567 if (index(item(2),"*").gt.0) then
571 ist2=isccortyp(rescode(1,item(2),0))
574 if (index(item(3),"*").gt.0) then
584 mask_tor(l,itypt1,itypt2,1)=masktemp
585 weitor(l,itypt1,itypt2,1)=weitemp
586 weitor_low(l,itypt1,itypt2,1)=weitemp_low
587 weitor_up(l,itypt1,itypt2,1)=weitemp_up
591 call card_concat(controlcard,.true.)
593 call card_concat(controlcard,.true.)
598 & write (iout,*) "Optimizable sidechain-torsional parameters:"
599 do itypt1=-nsccortyp,nsccortyp
600 do itypt2=-nsccortyp,nsccortyp
602 if (mask_tor(l,itypt1,itypt2,1).gt.0 .and. me.eq.Master)
603 & write (iout,'(4i5,3f10.5)') itypt1,itypt2,l,
604 & mask_tor(l,itypt1,itypt2,1),weitor(l,itypt1,itypt2,1),
605 & weitor_low(l,itypt1,itypt2,1),weitor_up(l,itypt1,itypt2,1)
611 c write (iout,*) "mod_fourier ",controlcard
613 mod_fourier(nloctyp)=index(controlcard,"MOD_FOURIER")
615 if (mod_fourier(nloctyp).gt.0) then
616 mod_fourier(nloctyp)=0
630 mask_eenew(ii,j,k,i)=0
638 call card_concat(controlcard,.true.)
639 do while ( index(controlcard,"END").eq.0 )
640 c write(iout,*) controlcard
641 call reads(controlcard,"TYPF",typf," ")
642 itypf=rescode(1,typf,0)
643 c write (iout,*) "TYPF ",typf," itypf",itypf
644 if (itypf.eq.0 .or. itypf.gt.ntyp) then
646 & write (*,*) "Error specifying fourier parms"
649 itypf=itype2loc(itypf)
650 c write (iout,*) "local type",itypf
652 if (index(controlcard,"B1_LOW").gt.0) then
654 call readi(controlcard,"IND",ind,0)
655 call readi(controlcard,"COEFF",ii,0)
656 if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
657 if (me.eq.Master) write (iout,*)
658 & "Error specifying B1, components undefined",ind,ii
661 mask_bnew1(ii,ind,itypf)=1
662 call reada(controlcard,"B1_LOW",bnew1_low(ii,ind,itypf),
664 call reada(controlcard,"B1_UP",bnew1_up(ii,ind,itypf),
666 mod_fourier(itypf)=mod_fourier(itypf)
667 & +mask_bnew1(ii,ind,itypf)
669 else if (index(controlcard,"B2_LOW").gt.0) then
671 mask_bnew2(ii,ind,itypf)=1
672 call readi(controlcard,"IND",ind,0)
673 call readi(controlcard,"COEFF",ii,0)
674 if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
675 if (me.eq.Master) write (iout,*)
676 & "Error specifying B2, components undefined",ind,ii
679 call reada(controlcard,"B2_LOW",bnew2_low(ii,ind,itypf),
681 call reada(controlcard,"B2_UP",bnew2_up(ii,ind,itypf),
683 mod_fourier(itypf)=mod_fourier(itypf)
684 & +mask_bnew2(ii,ind,itypf)
686 else if (index(controlcard,"C_LOW").gt.0) then
688 call readi(controlcard,"IND",ind,0)
689 call readi(controlcard,"COEFF",ii,0)
690 if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
691 if (me.eq.Master) write (iout,*)
692 & "Error specifying C, components undefined",ind,ii
695 mask_ccnew(ii,ind,itypf)=1
696 call reada(controlcard,"C_LOW",ccnew_low(ii,ind,itypf),.1d0)
697 call reada(controlcard,"C_UP",ccnew_up(ii,ind,itypf),0.0d0)
698 mod_fourier(itypf)=mod_fourier(itypf)
699 & +mask_ccnew(ii,ind,itypf)
701 else if (index(controlcard,"D_LOW").gt.0) then
703 call readi(controlcard,"IND",ind,0)
704 call readi(controlcard,"COEFF",ii,0)
705 if (ind.eq.0 .or. ii.eq.0 .or. ii.gt.3) then
706 if (me.eq.Master) write (iout,*)
707 & "Error specifying D, components undefined",ind,ii
710 mask_ddnew(ii,ind,itypf)=1
711 call reada(controlcard,"D_LOW",ddnew_low(ii,ind,itypf),
713 call reada(controlcard,"D_UP",ddnew_up(ii,ind,itypf),
715 mod_fourier(itypf)=mod_fourier(itypf)
716 & +mask_ddnew(ii,ind,itypf)
718 else if (index(controlcard,"E0_LOW").gt.0) then
720 call readi(controlcard,"COEFF",ii,0)
721 if (ii.eq.0 .or. ii.gt.3) then
722 if (me.eq.Master) write (iout,*)
723 & "Error specifying E0, components undefined",ii
726 mask_e0new(ii,itypf)=1
727 call reada(controlcard,"E0_LOW",e0new_low(ii,itypf),
729 call reada(controlcard,"E0_UP",e0new_up(ii,itypf),
731 mod_fourier(itypf)=mod_fourier(itypf)
732 & +mask_e0new(ii,itypf)
734 else if (index(controlcard,"E_LOW").gt.0) then
736 call readi(controlcard,"IND1",ind1,0)
737 call readi(controlcard,"IND2",ind2,0)
738 call readi(controlcard,"COEFF",ii,0)
739 if (ind1.eq.0 .or. ind2.eq.0 .or. ii.eq.0 .or. ii.gt.2) then
740 if (me.eq.Master) write (iout,*)
741 & "Error specifying E, components undefined",ind1,ind2,ii
744 mask_eenew(ii,ind1,ind2,itypf)=1
745 call reada(controlcard,"E_LOW",
746 & eenew_low(ii,ind1,ind2,itypf),0.1d0)
747 call reada(controlcard,"E_UP",
748 & eenew_up(ii,ind1,ind2,itypf),0.0d0)
749 mod_fourier(itypf)=mod_fourier(itypf)
750 & +mask_eenew(ii,ind1,ind2,itypf)
752 call card_concat(controlcard,.true.)
754 call card_concat(controlcard,.true.)
757 c write (iout,*) "itypf",itypf," mod_fourier",mod_fourier(itypf)
758 mod_fourier(nloctyp)=mod_fourier(nloctyp)+mod_fourier(itypf)
761 & write (iout,*) "mod_fourier all",mod_fourier(nloctyp)
763 call card_concat(controlcard,.true.)
764 do while ( index(controlcard,"END").eq.0 )
765 call readi(controlcard,"TYPF",typf," ")
767 if (itypf.eq.0 .or. itypf.gt.ntyp) then
769 & write (*,*) "Error specifying fourier parms"
772 itypf=itype2loc(itypf)
773 call readi(controlcard,"IND",ind,0)
774 call reada(controlcard,"B_LOW",b_low(ind,itypf),0.1d0)
775 call reada(controlcard,"B_UP",b_up(ind,itypf),0.0d0)
776 mask_fourier(ind,itypf)=1
777 mod_fourier(itypf)=mod_fourier(itypf)
778 & +mask_fourier(ind,itypf)
779 mod_fourier(nloctyp)=mod_fourier(nloctyp)
780 & +mask_fourier(ind,itypf)
781 call card_concat(controlcard,.true.)
783 call card_concat(controlcard,.true.)
793 mod_elec=index(controlcard,"MOD_ELEC").gt.0
796 call card_concat(controlcard,.true.)
797 do while ( index(controlcard,"END").eq.0 )
798 call readi(controlcard,"TYPE1",itype1,0)
799 call readi(controlcard,"TYPE2",itype2,0)
800 if (itype1.eq.0 .or. itype1.gt.2 .or. itype2.eq.0
801 & .or. itype2.gt.2) then
803 & write (iout,*) "Error specifying elec parms"
806 if (index(controlcard,"EPP").gt.0) then
808 mask_elec(itype1,itype2,1)=1
809 mask_elec(itype2,itype1,1)=1
810 call reada(controlcard,"EPP_LOW",epp_low(itype1,itype2),
812 epp_low(itype2,itype1)=epp_low(itype1,itype2)
813 call reada(controlcard,"EPP_UP",epp_up(itype1,itype2),
815 epp_up(itype2,itype1)=epp_up(itype1,itype2)
817 if (index(controlcard,"RPP").gt.0) then
819 mask_elec(itype1,itype2,2)=1
820 mask_elec(itype2,itype1,2)=1
821 call reada(controlcard,"RPP_LOW",rpp_low(itype1,itype2),
823 rpp_low(itype2,itype1)=rpp_low(itype1,itype2)
824 call reada(controlcard,"RPP_UP",rpp_up(itype1,itype2),
826 rpp_up(itype2,itype1)=rpp_up(itype1,itype2)
828 if (index(controlcard,"ELPP6").gt.0) then
830 mask_elec(itype1,itype2,3)=1
831 mask_elec(itype2,itype1,3)=1
832 call reada(controlcard,"ELPP6_LOW",
833 & elpp6_low(itype1,itype2),0.1d0)
834 elpp6_low(itype2,itype1)=elpp6_low(itype1,itype2)
835 call reada(controlcard,"ELPP6_UP",
836 & elpp6_up(itype1,itype2),0.0d0)
837 elpp6_up(itype2,itype1)=elpp6_up(itype1,itype2)
839 if (index(controlcard,"ELPP3").gt.0) then
841 mask_elec(itype1,itype2,4)=1
842 mask_elec(itype2,itype1,4)=1
843 call reada(controlcard,"ELPP3_LOW",
844 & elpp3_low(itype1,itype2),0.1d0)
845 elpp3_low(itype2,itype1)=elpp3_low(itype1,itype2)
846 call reada(controlcard,"ELPP3_UP",
847 & elpp3_up(itype1,itype2),0.0d0)
848 elpp3_up(itype2,itype1)=elpp3_up(itype1,itype2)
850 call card_concat(controlcard,.true.)
852 call card_concat(controlcard,.true.)
861 mod_scp=index(controlcard,"MOD_SCP").gt.0
864 call card_concat(controlcard,.true.)
865 do while (index(controlcard,"END").eq.0)
866 if (index(controlcard,"EPSCP").gt.0) then
868 call readi(controlcard,"ITYPSC",itypsc,0)
869 call readi(controlcard,"ITYPP",itypp,0)
870 if (itypsc.gt.20 .or. itypp.eq.0 .or. itypp.gt.2) then
872 & write (iout,*) "Error specifying scp parms"
875 mask_scp(itypsc,itypp,1)=1
876 call reada(controlcard,"EPSCP_LOW",
877 & epscp_low(itypsc,itypp),0.1d0)
878 call reada(controlcard,"EPSCP_UP",
879 & epscp_up(itypsc,itypp),0.0d0)
881 if (index(controlcard,"RSCP").gt.0) then
883 call readi(controlcard,"ITYPSC",itypsc,0)
884 call readi(controlcard,"ITYPP",itypp,0)
885 mask_scp(itypsc,itypp,2)=1
886 call readi(controlcard,"ITYPSC",itypsc,0)
887 call readi(controlcard,"ITYPP",itypp,0)
888 if (itypsc.gt.20 .or. itypp.eq.0 .or. itypp.gt.2) then
890 & write (iout,*) "Error specifying scp parms"
893 call reada(controlcard,"RSCP_LOW",
894 & rscp_low(itypsc,itypp),0.1d0)
895 c write(iout,*)itypsc,itypp,rscp_low(itypsc,itypp)
896 call reada(controlcard,"RSCP_UP",
897 & rscp_up(itypsc,itypp),0.0d0)
898 c write(iout,*)itypsc,itypp,rscp_up(itypsc,itypp)
900 call card_concat(controlcard,.true.)
903 c write (iout,*) "END ",controlcard
906 & mod_fourier(nloctyp).gt.0 .or. mod_elec .or. mod_scp
908 if (read_stat.lt.2. .and. mod_other_params) then
910 & write (iout,*)"Error - only weights and sidechain parameters",
911 & " can be optimized with READ_STAT=",read_stat
915 init_ene = init_ene .or. read_stat.eq.2
917 c write (*,*) "MOD_OTHER_PARAMS ",mod_other_params
918 c write (*,*) "MOD_SIDE ",mod_side," MOD_FOURIER",
919 c & mod_fourier(nloctyp),
920 c & " MOD_ELEC ",mod_elec," MOD_SCP ",mod_scp
921 init_ene=init_ene .or. mod_other_params
922 c write (iout,*) "init_ene",init_ene
927 c-----------------------------------------------------------------------------
928 subroutine print_general_data(*)
931 include "DIMENSIONS.ZSCOPT"
935 integer ierror,kolor,klucz
937 include "COMMON.WEIGHTS"
938 include "COMMON.NAMES"
939 include "COMMON.VMCPAR"
940 include "COMMON.TORSION"
941 include "COMMON.INTERACT"
942 include "COMMON.ENERGIES"
943 include "COMMON.MINPAR"
944 include "COMMON.IOUNITS"
945 include "COMMON.TIME1"
946 include "COMMON.PROTFILES"
947 include "COMMON.CHAIN"
948 include "COMMON.CLASSES"
949 include "COMMON.THERMAL"
950 include "COMMON.FFIELD"
951 include "COMMON.OPTIM"
952 include "COMMON.CONTROL"
953 include "COMMON.SCCOR"
954 character*800 controlcard
955 integer i,j,k,l,ii,n_ene_found,itypt,jtypt
956 integer ind,itype1,itype2,itypf,itypsc,itypp
957 integer ilen,lenpot,lenpre
959 character*3 liczba,liczba1
962 & write (iout,*) "Single-point target function evaluation"
963 else if (mode.eq.2) then
965 & write (iout,*) "Grid search of the parameter space"
966 else if (mode.eq.3) then
968 & write (iout,*) "Minimization of the target function"
971 & write (iout,*) "Wrong MODE",mode,", should be 1, 2, or 3"
975 if (me.eq.Master) write (iout,*) "RESCALE_MODE",rescale_mode
976 mod_other_params=index(controlcard,"OPTIMIZE_OTHER").gt.0
977 if (read_stat.eq.0 .and. mod_other_params) then
978 if (me.eq.Master) then
979 write (iout,*) "Error: only optimization of energy-term",
980 & " weights can be performed with READ_STAT=",read_stat
985 if (mode.eq.3 .and. me.eq.Master) then
986 write (iout,*) "Parameters of the SUMSL procedure:"
987 write (iout,'(a,t15,i5)') "MAXMIN",maxmin
988 write (iout,'(a,t15,i5)') "MAXFUN",maxfun
989 write (iout,'(a,t15,e15.7)') "TOLF",tolf
990 write (iout,'(a,t15,e15.7)') "RTOLF",rtolf
992 if (mod_other_params) then
993 if (me.eq.Master) write (iout,*)
994 &"Internal parameters of UNRES energy components will be optimized"
995 call card_concat(controlcard,.true.)
996 mod_side = (index(controlcard,"MOD_SIDE").gt.0)
997 if (mod_side .and. me.eq.Master) then
998 write (iout,*) "SC epsilons to be optimized:"
999 write (iout,*) "Single: eps(i,j)=eps(i,j)+(x(i)+x(j))/2"
1001 write (iout,*) restyp(ityp_ssc(i)),epss_low(i),epss_up(i)
1003 write (iout,*) "Pairs: eps(i,j)=eps(i,j)+x(i,j):"
1005 write (iout,*) restyp(ityp_psc(1,i)),
1006 & restyp(ityp_psc(2,i)),epsp_low(i),epsp_up(i)
1009 if (mod_sccor .and. me.eq.Master) then
1010 write (iout,*)"Torsional weights/coefficients to be optimized"
1011 write(iout,'(a)') "TYP1 TYP2 WEIGHT",
1012 & " LOWER BOUND UPPER BOUND"
1013 do itypt=-nsccortyp,nsccortyp
1014 do jtypt=-nsccortyp,nsccortyp
1016 if (mask_tor(l,itypt,jtypt,1).gt.0) then
1017 write(iout,'(3a4,3f10.5)')l,restyp(itypt),restyp(jtypt),
1018 & weitor(l,itypt,jtypt,1),weitor_low(l,itypt,jtypt,1),
1019 & weitor_up(l,itypt,jtypt,1)
1025 if (mod_fourier(nloctyp).gt.0 .and. me.eq.Master) then
1026 write (iout,*) "Fourier coefficients to be optimized"
1027 do itypf=0,nloctyp-1
1028 if (mod_fourier(itypf).gt.0) then
1029 write (iout,'(3a,i2)') "Type ",restyp(iloctyp(itypf)),
1030 & " number of coeffs",mod_fourier(itypf)
1031 write(iout,'(a)') "NAME COEFF LOWER BOUND UPPER BOUND"
1035 if (mask_bnew1(k,j,itypf).gt.0)
1036 & write (iout,'(2hB1,2i2,f10.5)') k,j,bnew1(k,j,itypf),
1037 & bnew1_low(k,j,itypf),bnew1_up(k,j,itypf)
1042 if (mask_bnew2(k,j,itypf).gt.0)
1043 & write (iout,'(2hB2,2i2,3f11.5)') k,j,bnew2(k,j,itypf),
1044 & bnew2_low(k,j,itypf),bnew2_up(k,j,itypf)
1049 if (mask_ccnew(k,j,itypf).gt.0)
1050 & write (iout,'(2hCC,2i2,3f11.5)') k,j,ccnew(k,j,itypf),
1051 & ccnew_low(k,j,itypf),ccnew_up(k,j,itypf)
1056 if (mask_ddnew(k,j,itypf).gt.0)
1057 & write (iout,'(2hDD,2i2,3f11.5)') k,j,ddnew(k,j,itypf),
1058 & ddnew_low(k,j,itypf),ddnew_up(k,j,itypf)
1062 if (mask_e0new(j,itypf).gt.0)
1063 & write (iout,'(2hE0,i2,3f11.5)') j,e0new(j,itypf),
1064 & e0new_low(j,itypf),e0new_up(j,itypf)
1069 if (mask_eenew(l,k,j,itypf).gt.0)
1070 & write (iout,'(2hEE,3i2,3f11.5)')
1071 & l,k,j,eenew(l,k,j,itypf),eenew_low(l,k,j,itypf),
1072 & eenew_up(l,k,j,itypf)
1078 if (mask_fourier(i,itypf).gt.0) then
1079 write (iout,'(1hB,i2,3f11.5)')
1080 & i,b(i,itypf),b_low(i,itypf),b_up(i,itypf)
1087 if (mod_elec .and. Me.eq.Master) then
1089 write (iout,*) "Electrostatic parameters to be optimized"
1092 if (mask_elec(itype1,itype2,1).gt.0)
1093 & write (iout,'(2i3," EPP",f8.3," EPP_LOW",f8.3,
1095 & itype1,itype2,epp(itype1,itype2),epp_low(itype1,itype2),
1096 & epp_up(itype1,itype2)
1097 if (mask_elec(itype1,itype2,2).gt.0)
1098 & write (iout,'(2i3," RPP",f8.3," RPP_LOW",f8.3,
1100 & itype1,itype2,rpp(itype1,itype2),rpp_low(itype1,itype2),
1101 & rpp_up(itype1,itype2)
1102 if (mask_elec(itype1,itype2,3).gt.0)
1103 & write (iout,'(2i3," ELPP6",f8.3," ELPP6_LOW",f8.3,
1104 & " ELPP6_UP",f8.3)')
1105 & itype1,itype2,elpp6(itype1,itype2),
1106 & elpp6_low(itype1,itype2),elpp6_up(itype1,itype2)
1107 if (mask_elec(itype1,itype2,4).gt.0)
1108 & write (iout,'(2i3," ELPP3",f8.3," ELPP3_LOW",f8.3,
1109 & " ELPP3_UP",f8.3)')
1110 & itype1,itype2,elpp3(itype1,itype2),
1111 & elpp3_low(itype1,itype2),elpp3_up(itype1,itype2)
1115 if (mod_scp .and. me.eq.Master) then
1118 write (iout,*) "SCP parameters to be optimized:"
1119 if (mask_scp(0,1,1)+mask_scp(0,2,1)+mask_scp(0,1,2)+
1120 & mask_scp(0,2,2) .gt. 0) then
1122 & "SCP parameters are assumed to depend on peptide group type only"
1124 if (mask_scp(0,j,1).gt.0)
1125 & write (iout,'(i3," EPSCP",f8.3," EPSCP_LOW",f8.3,
1126 & " EPSCP_UP",f8.3)') j,eps_scp(1,j),epscp_low(0,j),
1128 if (mask_scp(0,j,2).gt.0)
1129 & write (iout,'(i3," RSCP",f8.3," RSCP_LOW",f8.3,
1130 & " RSCP_UP",f8.3)') j,rscp(1,j),rscp_low(0,j),
1136 if (mask_scp(i,j,1).gt.0)
1137 & write (iout,'(2i3," EPSCP",f8.3," EPSCP_LOW",f8.3,
1138 & " EPSCP_UP",f8.3)') i,j,eps_scp(i,j),epscp_low(i,j),
1140 if (mask_scp(i,j,2).gt.0)
1141 & write (iout,'(2i3," RSCP",f8.3," RSCP_LOW",f8.3,
1142 & " RSCP_UP",f8.3)') i,j,rscp(i,j),rscp_low(i,j),
1150 & write (iout,*) "MOD_OTHER_PARAMS ",mod_other_params
1152 & write (iout,*) "MOD_SIDE ",mod_side," MOD_FOURIER",
1153 & mod_fourier(nloctyp),
1154 & " MOD_ELEC ",mod_elec," MOD_SCP ",mod_scp
1155 init_ene=init_ene .or. mod_other_params
1156 if (me.eq.Master) write (iout,*) "init_ene",init_ene
1157 if (me.eq.Master) write (iout,*)
1161 c-----------------------------------------------------------------------------
1162 subroutine read_protein_data(*)
1164 include "DIMENSIONS"
1165 include "DIMENSIONS.ZSCOPT"
1168 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1169 include "COMMON.MPI"
1171 include "COMMON.CONTROL"
1172 include "COMMON.CHAIN"
1173 include "COMMON.CLASSES"
1174 include "COMMON.COMPAR"
1175 include "COMMON.ENERGIES"
1176 include "COMMON.IOUNITS"
1177 include "COMMON.PROTFILES"
1178 include "COMMON.PROTNAME"
1179 include "COMMON.VMCPAR"
1180 include "COMMON.OPTIM"
1181 include "COMMON.WEIGHTS"
1182 include "COMMON.XBOUND"
1183 include "COMMON.NAMES"
1184 include "COMMON.ALLPROT"
1185 include "COMMON.THERMAL"
1186 include "COMMON.TIME1"
1187 include "COMMON.WEIGHTDER"
1188 character*64 nazwa,key
1189 character*16000 controlcard
1190 character*16000 all_protfiles
1191 integer i,j,k,kk,l,iprot,ii,if,ib,ibatch,nn,nn1,iww,maskcheck,
1192 & il,inat,ilevel,weightread,jfrag,jclass,nfragm,iparm
1193 integer nrec,nlines,iscor
1194 double precision energ,wangnorm(maxT),wq(maxT),wrms(maxT),
1195 & wrgy(maxT),wsign(maxT),wangnat(maxT),wqnat(maxT),wrmsnat(maxT),
1196 & wrgynat(maxT),wcorangnorm(maxT),wcorq(maxT),
1197 & wcorrms(maxT),wcorrgy(maxT),wcorsign(maxT),wcorangnat(maxT),
1198 & wcorqnat(maxT),wcorrmsnat(maxT),wcorrgynat(maxT),
1199 & angnormlow(maxT),qlow(maxT),rmslow(maxT),
1200 & rgylow(maxT),signlow(maxT),angnatlow(maxT),
1201 & qnatlow(maxT),rmsnatlow(maxT),rgynatlow(maxT),
1202 & angcorlow(maxT),qcorlow(maxT),rmscorlow(maxT),rgycorlow(maxT),
1203 & signcorlow(maxT),angcornatlow(maxT),qcornatlow(maxT),
1204 & rmscornatlow(maxT),rgycornatlow(maxT),signcornatlow(maxT),
1205 & angnormup(maxT),qup(maxT),rmsup(maxT),rgyup(maxT),signup(maxT),
1206 & angnatup(maxT),qnatup(maxT),rmsnatup(maxT),rgynatup(maxT),
1207 & angcorup(maxT),qcorup(maxT),rmscorup(maxT),rgycorup(maxT),
1209 & angcornatup(maxT),qcornatup(maxT),rmscornatup(maxT),
1210 & rgycornatup(maxT),signcornatup(MaxT)
1213 double precision ebia(maxprot),rjunk
1215 character*64 zeros /
1216 &'0000000000000000000000000000000000000000000000000000000000000000'
1219 c print *,"Processor",me," calls read_protein_data"
1221 C Read seventh record: general data of proteins used for calibration
1222 call card_concat(controlcard,.true.)
1223 c write(2, *)controlcard
1224 call readi(controlcard,"NPROT",nprot,0)
1225 pdbref=index(controlcard,"PDBREF").gt.0
1226 print_refstr=me.eq.Master.and.
1227 & index(controlcard,"PRINT_REFSTR").gt.0
1228 if (nprot.eq.0) then
1229 if (me.eq.Master) then
1230 write(iout,*) "Error: at least one protein must be specified."
1237 read (inp,'(a)') protname(iprot)
1238 if (me.eq.Master) then
1240 write (iout,*) "Reading data of protein",iprot," named ",
1243 call card_concat(controlcard,.true.)
1244 call reada(controlcard,"ENECUT_MIN",enecut_min(iprot),15.0d0)
1245 call reada(controlcard,"ENECUT_MAX",enecut_max(iprot),100.0d0)
1246 call reada(controlcard,"ENECUT",enecut(iprot),enecut_min(iprot))
1247 if (enecut(iprot).lt.enecut_min(iprot))
1248 & enecut(iprot)=enecut_min(iprot)
1249 if (enecut_max(iprot).le.enecut_min(iprot))
1250 & enecut_max(iprot)=2*enecut_min(iprot)
1252 & write (iout,'(3(a,f9.1))') "ENECUT",enecut(iprot)," ENECUT_MIN",
1253 & enecut_min(iprot)," ENECUT_MAX",enecut_max(iprot)
1254 call readi(controlcard,"HEFAC",hefac(iprot),50)
1255 call readi(controlcard,"HTFAC",htfac(iprot),50)
1256 call readi(controlcard,"HELOW",hemax_low(iprot),20)
1257 call readi(controlcard,"HTLOW",htmax_low(iprot),20)
1258 if (me.eq.Master) write (iout,*) "iprot",iprot,
1259 & " hefac",hefac(iprot)," helow",hemax_low(iprot),
1260 & " htfac",htfac(iprot)," htlow",htmax_low(iprot)
1261 c 7/27/2013 AL Read maximum likelihood data
1262 call card_concat(controlcard,.true.)
1263 call readi(controlcard,"NBETA_L",nbeta(1,iprot),0)
1264 if (me.eq.Master) write (iout,*) "NBETA_L",nbeta(1,iprot)
1265 caonly(iprot)=index(controlcard,"CAONLY").gt.0
1266 sconly(iprot)=index(controlcard,"SCONLY").gt.0
1267 rmscomp(iprot)=index(controlcard,"RMSCOMP").gt.0
1268 call reada(controlcard,"SIGMA",sigma2(iprot),4.0d0)
1270 & write (iout,*) "RMSCOMP",rmscomp," SIGMA",sigma2(iprot),
1271 & " CAONLY ",caonly(iprot)," SCONLY",sconly(iprot)
1272 do ib=1,nbeta(1,iprot)
1273 read(inp,*)betaT(ib,1,iprot),weilik(ib,iprot),
1276 & write (iout,*) i,betaT(ib,1,iprot),weilik(ib,iprot),
1279 c 7/27/2013 AL Read heat-capacity data
1280 call card_concat(controlcard,.true.)
1281 call readi(controlcard,"NBETA_CV",nbeta(2,iprot),0)
1282 call reada(controlcard,"WCV",wcv(iprot),1.0d0)
1283 call reada(controlcard,"BASE",heatbase(iprot),0.0d0)
1285 & write (iout,*) "NBETA_CV",nbeta(2,iprot)," WCV",wcv(iprot)
1286 do ib=1,nbeta(2,iprot)
1287 read(inp,*) betaT(ib,2,iprot),target_cv(ib,iprot),
1289 weiCv(ib,iprot)=weiCv(ib,iprot)*wcv(iprot)
1291 & write (iout,*) betaT(ib,2,iprot),target_cv(ib,iprot),
1294 if (me.eq.Master) then
1295 write (iout,*) "Experimental heat capacity curve"
1296 do ib=1,nbeta(2,iprot)
1297 write (iout,'(f6.2,2f10.5)') betaT(ib,2,iprot),
1298 & target_cv(ib,iprot),weiCv(ib,iprot)
1300 write (iout,'(a,f10.5)') "Baseline",heatbase(iprot)
1303 c Conformation-dependent averages
1304 call card_concat(controlcard,.true.)
1305 call readi(controlcard,"NATLIKE",natlike(iprot),0)
1306 do i=1,natlike(iprot)
1307 call card_concat(controlcard,.true.)
1308 call reada(controlcard,"WNAT",wnat(i,iprot),1.0d0)
1309 call readi(controlcard,"NUMNAT",numnat(i,iprot),1)
1310 call readi(controlcard,"NATDIM",natdim(i,iprot),1)
1311 do ib=1,nbeta(i+2,iprot)
1312 read(inp,*)betaT(ib,i+2,iprot),(weinat(k,ib,i,iprot),
1313 & nuexp(k,ib,i,iprot),k=1,natdim(i,iprot))
1316 do i=1,natlike(iprot)+2
1317 do j=1,nbeta(i,iprot)
1318 betaT(j,i,iprot)=1.0d0/(Rgas*betaT(j,i,iprot))
1320 & write (2,*) "R i",i," j",j," beta",betaT(j,i,iprot)
1325 C Read names of files with the data for protein IPROT
1326 call card_concat(controlcard,.false.)
1328 if (iparm.eq.myparm) then
1329 call split_string(controlcard,protfiles(1,iprot),
1330 & maxfile_prot,nfile_prot(iprot))
1332 write(iout,*)"iprot",iprot," nfile",nfile_prot(iprot)
1334 & (protfiles(i,iprot),i=1,nfile_prot(iprot))
1340 c Read molecular information of the protein
1341 call molread_zs(iprot)
1342 c 3/31/04 AL Read the reference structures of protein iprot
1343 c print *,"Calling read_ref_structure"
1344 call read_ref_structure(iprot,*11)
1346 write (iout,*) "Protein",iprot," left READ_REF_STRUCTURE"
1352 c-------------------------------------------------------------------------------
1353 subroutine read_database(*)
1355 include "DIMENSIONS"
1356 include "DIMENSIONS.ZSCOPT"
1359 integer IERROR,ERRCODE,STATUS(MPI_STATUS_SIZE)
1360 include "COMMON.MPI"
1362 include "COMMON.CHAIN"
1363 include "COMMON.INTERACT"
1364 include "COMMON.CLASSES"
1365 include "COMMON.ENERGIES"
1366 include "COMMON.IOUNITS"
1367 include "COMMON.PROTFILES"
1368 include "COMMON.PROTNAME"
1369 include "COMMON.VMCPAR"
1370 include "COMMON.WEIGHTS"
1371 include "COMMON.NAMES"
1372 include "COMMON.ALLPROT"
1373 include "COMMON.WEIGHTDER"
1374 include "COMMON.VAR"
1375 include "COMMON.SBRIDGE"
1376 include "COMMON.GEO"
1377 include "COMMON.COMPAR"
1378 include "COMMON.OPTIM"
1380 character*16000 controlcard
1381 character*16000 all_protfiles
1382 character*3 liczba,liczba1
1383 integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,if,ib,ibatch,
1385 integer ixdrf,iret,itmp
1386 integer nrec,nlines,iscor
1387 double precision energ,t_acq,tcpu
1390 double precision rjunk
1391 integer ntot_all(0:maxprot,0:maxprocs-1)
1393 double precision energia(0:max_ene),etot
1394 real*4 csingle(3,maxres2),reini,refree,rmsdev,prec
1395 integer Previous,Next
1396 c print *,"Processor",me," calls read_protein_data"
1398 if (me.eq.master) then
1399 Previous=MPI_PROC_NULL
1403 if (me.eq.nprocs-1) then
1409 c Set the scratchfile names
1410 write (liczba,'(bz,i3.3)') me
1411 write (liczba1,'(bz,i3.3)') myparm
1413 c 1/27/05 AL Change stored coordinates to single precision and don't store
1414 c energy components in the binary databases.
1415 lenrec(iprot)=12*(nres_zs(iprot)+nct_zs(iprot)-nnt_zs(iprot)+1)
1416 & +4*(2*nss_zs(1,iprot)+1)+8*natlike(iprot)*maxdimnat+16
1417 c 4/13/04 AL Add space for similarity measure
1418 lenrec_ene(iprot) = (2*nntyp+3*n_ene+2)*8
1419 & +8*natlike(iprot)*maxdimnat
1422 write (iout,*) "Protein i"," lenrec",lenrec(iprot)
1423 write (iout,*) "lenrec_ene",lenrec_ene(iprot)
1426 bprotfiles(iprot) = scratchdir(:ilen(scratchdir))//
1427 & "/"//protname(iprot)(:ilen(protname(iprot)))//
1428 & "_par"//liczba1//"_"//liczba//".xbin"
1429 benefiles(iprot) = scratchdir(:ilen(scratchdir))//
1430 & "/"//protname(iprot)(:ilen(protname(iprot)))//
1431 & "_par"//liczba1//"_"//liczba//".enebin"
1432 c write (iout,*) "scratchfile ",
1433 c & bprotfiles(iprot)(:ilen(bprotfiles(iprot)))
1442 call restore_molinfo(iprot)
1443 open (ientout,file=bprotfiles(iprot),status="unknown",
1444 & form="unformatted",access="direct",recl=lenrec(iprot))
1445 open (istat,file=benefiles(iprot),status="unknown",
1446 & form="unformatted",access="direct",recl=lenrec_ene(iprot))
1447 c Read conformations from binary DA files (one per batch) and write them to
1448 c a binary DA scratchfile.
1451 write (liczba,'(bz,i3.3)') me
1453 IF (ME.EQ.MASTER) THEN
1454 c Only the master reads the database; it'll send it to the other procs
1460 do if=1,nfile_prot(iprot)
1461 nazwa=protfiles(if,iprot)
1462 & (:ilen(protfiles(if,iprot)))//".cx"
1464 write (iout,*) "Opening file ",nazwa(:ilen(nazwa))
1466 #if (defined(AIX) && !defined(JUBL))
1467 call xdrfopen_(ixdrf,nazwa, "r", iret)
1469 call xdrfopen(ixdrf,nazwa, "r", iret)
1471 if (iret.eq.0) goto 1111
1475 #if (defined(AIX) && !defined(JUBL))
1476 call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret)
1477 if (iret.eq.0) goto 102
1478 call xdrfint_(ixdrf, nss, iret)
1479 if (iret.eq.0) goto 102
1481 call xdrfint_(ixdrf, ihpb(j), iret)
1482 if (iret.eq.0) goto 102
1483 call xdrfint_(ixdrf, jhpb(j), iret)
1484 if (iret.eq.0) goto 102
1486 call xdrffloat_(ixdrf,reini,iret)
1487 if (iret.eq.0) goto 102
1488 call xdrffloat_(ixdrf,refree,iret)
1489 if (iret.eq.0) goto 102
1490 call xdrfint_(ixdrf,natlik,iret)
1491 if (iret.eq.0) goto 102
1493 call xdrfint(ixdrf,natliktemp(j),iret)
1494 if (iret.eq.0) goto 102
1495 do k=1,natliktemp(j)
1496 call xdrffloat(ixdrf,nutemp(k,j),iret)
1497 if (iret.eq.0) goto 102
1501 c write (0,*) "me",me," iprot",iprot," i",i
1502 call xdrf3dfcoord(ixdrf, csingle, itmp, prec, iret)
1503 if (iret.eq.0) goto 102
1504 call xdrfint(ixdrf, nss, iret)
1505 if (iret.eq.0) goto 102
1507 call xdrfint(ixdrf, ihpb(k), iret)
1508 if (iret.eq.0) goto 102
1509 call xdrfint(ixdrf, jhpb(k), iret)
1510 if (iret.eq.0) goto 102
1512 call xdrffloat(ixdrf,reini,iret)
1513 if (iret.eq.0) goto 102
1514 call xdrffloat(ixdrf,refree,iret)
1515 if (iret.eq.0) goto 102
1517 call xdrfint(ixdrf,natlik,iret)
1518 if (iret.eq.0) goto 102
1520 call xdrfint(ixdrf,natliktemp(j),iret)
1521 if (iret.eq.0) goto 102
1522 do k=1,natliktemp(j)
1523 call xdrffloat(ixdrf,nutemp(k,j),iret)
1524 if (iret.eq.0) goto 102
1529 call xdrffloat(ixdrf,rmsdev,iret)
1530 if (iret.eq.0) goto 102
1531 c write (2,*) "rmsdev",rmsdev
1532 call xdrfint(ixdrf,iscor,iret)
1533 if (iret.eq.0) goto 102
1534 c write (2,*) "iscor",iscor
1537 eini(jj+1,iprot)=reini
1538 entfac(jj+1,iprot)=refree
1546 c(l,nres+k)=csingle(l,nres+k-nnt+1)
1550 write (iout,'(5hREAD ,i5,2f15.4)')
1551 & jj+1,eini(jj+1,iprot),entfac(jj+1,iprot)
1552 write (iout,*) "natlik",natlik
1554 write (iout,*) "natliketemp",natliktemp(l)
1555 write(iout,*) (nutemp(k,l),k=1,natliktemp(l))
1557 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1558 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1561 call add_new_cconf(jjj,jj,jj_old,icount,iprot,
1566 & write (iout,*) "Protein ",protname(iprot),
1567 & i-1," conformations read from DA file ",
1568 & nazwa(:ilen(nazwa))
1570 & write (iout,*) jj," conformations read so far"
1571 #if (defined(AIX) && !defined(JUBL))
1572 call xdrfclose_(ixdrf,nazwa,iret)
1574 call xdrfclose(ixdrf,nazwa,iret)
1576 c print *,"file ",nazwa(:ilen(nazwa))," closed"
1580 write (iout,*) "jj_old",jj_old," jj",jj
1582 call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
1583 if (icount.gt.0) call MPI_Send(0,1,MPI_INTEGER,Next,570,
1587 call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
1589 t_acq = tcpu() - t_acq
1591 & write (iout,*) "Processor",me," protein",iprot,
1592 & " batch",ibatch," time for conformation read/send",t_acq
1595 c A worker gets the confs from the master and sends them to its neighbor
1597 call receive_and_pass_cconf(icount,jj_old,jj,iprot,
1599 t_acq = tcpu() - t_acq
1600 c write (iout,*) "Processor",me," protein",iprot,
1601 c & " batch",ibatch,
1602 c & " time for conformation receive/send",t_acq
1606 if (me.eq.Master) write (iout,*) "Protein",
1607 & protname(iprot)(:ilen(protname(iprot))),", ",ntot(iprot),
1608 & " conformatons read ",ntot(iprot)," conformations written to ",
1609 & bprotfiles(iprot)(:ilen(bprotfiles(iprot)))
1610 ntot(0) = ntot(0)+ntot(iprot)
1616 & write(iout,*)"A total of",ntot(0)," conformations read."
1618 c Check if everyone has the same number of conformations
1619 call MPI_Allgather(ntot(0),maxprot+1,MPI_INTEGER,
1620 & ntot_all(0,0),maxprot+1,MPI_INTEGER,ALL_Comm,IERROR)
1625 if (ntot(j).ne.ntot_all(j,i)) then
1627 & write (iout,*) "Number of conformations at processor",i,
1628 & " for protein",j," differs from that at processor",me,
1629 & ntot(j),ntot_all(j,i)
1636 c----------- Temporary; reading probs from external file
1637 open (88,file='1LE1_wham_last_2.rms',status='old')
1639 read (88,*) ii,weirms(i,1)
1642 if (me.eq.Master) write (iout,*) "i",i," weirms",weirms(i,1)
1645 call MPI_Bcast(weirms(1,1), ntot(1), MPI_Double_Precision,
1646 & Master, ALL_COMM, IERROR)
1647 c----------- end temportary stuff
1650 if (me.eq.Master) then
1652 write (iout,*) "Number of conformations read by processors"
1653 write (iout,'(10x,7a10)') (protname(i),i=0,nprot)
1656 write (iout,'(8i10)') i,(ntot_all(j,i),j=0,nprot)
1658 write (iout,*) "Calculation terminated."
1665 1111 if (me.eq.Master)
1666 &write(iout,*) "Error opening coordinate file ",nazwa(:ilen(nazwa))
1669 call MPI_Abort(ALL_COMM,IERROR,ERRCODE)
1673 c------------------------------------------------------------------------------
1674 subroutine add_new_cconf(jjj,jj,jj_old,icount,iprot,Next)
1676 include "DIMENSIONS"
1677 include "DIMENSIONS.ZSCOPT"
1678 include "COMMON.CHAIN"
1679 include "COMMON.INTERACT"
1680 include "COMMON.LOCAL"
1681 include "COMMON.CLASSES"
1682 include "COMMON.ENERGIES"
1683 include "COMMON.IOUNITS"
1684 include "COMMON.PROTFILES"
1685 include "COMMON.PROTNAME"
1686 include "COMMON.VMCPAR"
1687 include "COMMON.WEIGHTS"
1688 include "COMMON.NAMES"
1689 include "COMMON.ALLPROT"
1690 include "COMMON.WEIGHTDER"
1691 include "COMMON.VAR"
1692 include "COMMON.SBRIDGE"
1693 include "COMMON.GEO"
1694 include "COMMON.COMPAR"
1698 integer i,j,jj,jjj,jj_old,icount,k,kk,l,iprot,ii,ib,ibatch,
1699 & nn,nn1,inan,Next,itj
1700 double precision etot,energia(0:max_ene)
1702 c if (entfac(jj+1,iprot).gt.-10.0d0
1703 c & .or. entfac(jj+1,iprot).lt.-150.0d0) then
1704 c write (iout,*) "Entropy factor out of range for conformation",
1705 c & jjj,entfac(jj+1,iprot),", conformation skipped."
1708 call int_from_cart1(.false.)
1710 if (vbld(j).lt.2.0d0 .or. vbld(j).gt.5.0d0) then
1711 write (iout,*) "nnt",nnt," nct",nct
1712 write (iout,*) "Bad CA-CA bond length",j," ",vbld(j)
1713 write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
1714 write (iout,*) "The Cartesian geometry is:"
1715 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1716 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1717 write (iout,*) "The internal geometry is:"
1718 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1719 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1720 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1721 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1722 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1723 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1725 & "This conformation WILL NOT be added to the database."
1731 if (itype(j).ne.10 .and. (vbld(nres+j)-dsc(itj)).gt.2.0d0) then
1732 write (iout,*) "nnt",nnt," nct",nct
1733 write (iout,*) "Bad CA-SC bond length",ii," ",vbld(nres+j)
1734 write (iout,*) "The Cartesian geometry is:"
1735 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1736 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1737 write (iout,*) "The internal geometry is:"
1738 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1739 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1740 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1741 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1742 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1743 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1745 & "This conformation WILL NOT be added to the database."
1750 if (theta(j).le.0.0d0) then
1752 & "Zero theta angle(s) in conformation",ii
1753 write (iout,*) "The Cartesian geometry is:"
1754 write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres)
1755 write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct)
1756 write (iout,*) "The internal geometry is:"
1757 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1758 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1759 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1760 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1761 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1762 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1764 & "This conformation WILL NOT be added to the database."
1767 if (theta(j).gt.179.97*deg2rad) theta(j)=179.97*deg2rad
1769 if (.not. init_ene) then
1772 etot=etot+ww(j)*enetb(icount+1,j,iprot)
1778 if (isnan(etot).ne.0) inan=1
1780 if (isnan(real(etot))) inan=1
1784 idumm=proc_proc(etot,inan)
1786 call proc_proc(etot,inan)
1793 write (iout,*) "NaNs detected in some of the energy",
1794 & " components for protein",iprot," batch",ibatch,
1795 & " conformation",jjj
1796 write (iout,*) "etot",etot
1797 write (iout,*) "The Cartesian geometry is:"
1798 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
1799 write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
1800 write (iout,*) "The internal geometry is:"
1801 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1802 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1803 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1804 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1805 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1806 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1807 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1808 write (iout,*) "The components of the energy are:"
1811 energia(k)=enetb(jj+1,k,iprot)
1813 call enerprint(energia(0))
1815 & "This conformation WILL NOT be added to the database."
1820 write (iout,'(e15.5,16i5)') entfac(icount+1,iprot),
1821 & iscore(icount+1,0,iprot),ibatch
1822 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
1823 write (iout,'(8f10.5)') ((c(j,i+nres),j=1,3),i=nnt,nct)
1824 write (iout,'(8f10.4)') (vbld(k),k=nnt+1,nct)
1825 write (iout,'(8f10.4)') (vbld(k),k=nres+nnt,nres+nct)
1826 write (iout,'(8f10.4)') (rad2deg*theta(k),k=3,nres)
1827 write (iout,'(8f10.4)') (rad2deg*phi(k),k=4,nres)
1828 write (iout,'(8f10.4)') (vbld(k+nres),k=nnt,nct)
1829 write (iout,'(8f10.4)') (rad2deg*alph(k),k=2,nres-1)
1830 write (iout,'(8f10.4)') (rad2deg*omeg(k),k=2,nres-1)
1831 write (iout,'(16i5)') nss,(ihpb(k),jhpb(k),k=1,nss)
1832 write (iout,'(5e15.5)') (enetb(icount+1,j,iprot),j=1,n_ene)
1833 c write (iout,'(2e15.5)') (eneps(1,j,icount+1,iprot),
1834 c & eneps(2,j,icount+1,iprot),j=1,nntyp)
1836 c write (iout,*) "First nu in readrtms"
1839 list_conf(jj,iprot)=jj
1840 call store_cconf_from_file(jj,icount,iprot)
1841 do j=1,natlike(iprot)
1843 if (k.eq.numnat(j,iprot)) then
1844 do l=1,natdim(j,iprot)
1845 nu(l,k,j,iprot)=nutemp(l,k)
1851 c write (iout,*) "End First nu in readrtms"
1853 if (icount.eq.maxstr_proc) then
1855 write (iout,* ) "jj_old",jj_old," jj",jj
1856 write (iout,*) "Calling write_and_send_cconf"
1859 call write_and_send_cconf(icount,jj_old,jj,iprot,Next)
1862 write (iout,*) "Exited write_and_send_cconf"
1870 c------------------------------------------------------------------------------
1871 subroutine store_cconf_from_file(jj,icount,iprot)
1873 include "DIMENSIONS"
1874 include "DIMENSIONS.ZSCOPT"
1875 include "COMMON.CHAIN"
1876 include "COMMON.SBRIDGE"
1877 include "COMMON.INTERACT"
1878 include "COMMON.IOUNITS"
1879 include "COMMON.CLASSES"
1880 include "COMMON.ALLPROT"
1881 include "COMMON.VAR"
1882 integer i,j,jj,icount,ibatch,iprot
1883 c Store the conformation that has been read in
1886 c_zs(j,i,icount,iprot)=c(j,i)
1889 nss_zs(icount,iprot)=nss
1891 ihpb_zs(i,icount,iprot)=ihpb(i)
1892 jhpb_zs(i,icount,iprot)=jhpb(i)
1896 c------------------------------------------------------------------------------
1897 subroutine write_and_send_cconf(icount,jj_old,jj,iprot,Next)
1899 include "DIMENSIONS"
1900 include "DIMENSIONS.ZSCOPT"
1904 include "COMMON.MPI"
1906 include "COMMON.CHAIN"
1907 include "COMMON.SBRIDGE"
1908 include "COMMON.INTERACT"
1909 include "COMMON.IOUNITS"
1910 include "COMMON.CLASSES"
1911 include "COMMON.VAR"
1912 include "COMMON.ALLPROT"
1913 include "COMMON.ENERGIES"
1914 include "COMMON.WEIGHTDER"
1915 include "COMMON.OPTIM"
1916 include "COMMON.COMPAR"
1917 integer icount,jj_old,jj,Next,iprot
1918 c Write the structures to a scratch file
1920 c Master sends the portion of conformations that have been read in to the neighbor
1922 write (iout,*) "Processor",me," entered WRITE_AND_SEND_CONF"
1925 call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,IERROR)
1927 write (iout,*) "Processor",me," Next",next," sent icount=",icount
1928 write (iout,*) "Processor",me," jj_old",jj_old," jj",jj
1931 if (icount.gt.0) then
1932 call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER,
1933 & Next,571,WHAM_COMM,IERROR)
1934 call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
1935 & Next,572,WHAM_COMM,IERROR)
1936 if (.not. init_ene) call MPI_Send(enetb(jj_old,1,iprot),
1938 & MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,IERROR)
1939 call MPI_Send(nu(1,1,jj_old,iprot),
1940 & maxdimnat*natlike(iprot)*icount,
1941 & MPI_DOUBLE_PRECISION,
1942 & Next,577,WHAM_COMM,IERROR)
1943 call MPI_Send(eini(jj_old,iprot),icount,
1944 & MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR)
1945 call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
1946 & Next,579,WHAM_COMM,IERROR)
1947 call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2,
1948 & MPI_REAL,Next,580,WHAM_COMM,IERROR)
1949 if (.not. init_ene) call MPI_Send(eneps(1,1,1,iprot),
1951 & MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR)
1954 write (2,*) "Calling dawrite_ccoords"
1956 call dawrite_ccoords(iprot,jj_old,jj,ientout)
1957 call dawrite_ene(iprot,jj_old,jj,istat)
1960 c------------------------------------------------------------------------------
1962 subroutine receive_and_pass_cconf(icount,jj_old,jj,iprot,Previous,
1965 include "DIMENSIONS"
1966 include "DIMENSIONS.ZSCOPT"
1968 integer IERROR,STATUS(MPI_STATUS_SIZE)
1969 include "COMMON.MPI"
1970 include "COMMON.CHAIN"
1971 include "COMMON.SBRIDGE"
1972 include "COMMON.INTERACT"
1973 include "COMMON.IOUNITS"
1974 include "COMMON.CLASSES"
1975 include "COMMON.ALLPROT"
1976 include "COMMON.ENERGIES"
1977 include "COMMON.VAR"
1978 include "COMMON.GEO"
1979 include "COMMON.WEIGHTDER"
1980 include "COMMON.COMPAR"
1981 include "COMMON.OPTIM"
1982 integer i,j,k,icount,jj_old,jj,iprot,Previous,Next
1985 write (iout,*) "Processor",me," entered RECEIVE_AND_PASS_CONF"
1988 do while (icount.gt.0)
1989 c call MPI_Probe(Previous,570,WHAM_COMM,STATUS,IERROR)
1990 call MPI_Recv(icount,1,MPI_INTEGER,Previous,570,WHAM_COMM,
1993 write (iout,*)"Processor",me," previous",previous," icount",icount
1996 call MPI_Send(icount,1,MPI_INTEGER,Next,570,WHAM_COMM,
1999 write (iout,*) "Processor",me," icount",icount
2002 if (icount.eq.0) return
2003 call MPI_Recv(nss_zs(1,iprot),icount,MPI_INTEGER,
2004 & Previous,571,WHAM_COMM,STATUS,IERROR)
2005 call MPI_Send(nss_zs(1,iprot),icount,MPI_INTEGER,
2006 & Next,571,WHAM_COMM,IERROR)
2007 call MPI_Recv(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2008 & Previous,572,WHAM_COMM,STATUS,IERROR)
2009 call MPI_Send(ihpb_zs(1,1,iprot),icount,MPI_INTEGER,
2010 & Next,572,WHAM_COMM,IERROR)
2011 if (.not. init_ene) then
2012 call MPI_Recv(enetb(jj_old,1,iprot),maxstr*n_ene,
2013 & MPI_DOUBLE_PRECISION,Previous,576,WHAM_COMM,STATUS,IERROR)
2014 call MPI_Send(enetb(jj_old,1,iprot),maxstr*n_ene,
2015 & MPI_DOUBLE_PRECISION,Next,576,WHAM_COMM,STATUS,IERROR)
2017 call MPI_Recv(nu(1,1,jj_old,iprot),
2018 & maxdimnat*natlike(iprot)*icount,
2019 & MPI_DOUBLE_PRECISION,
2020 & Previous,577,WHAM_COMM,STATUS,IERROR)
2021 call MPI_Send(nu(1,1,jj_old,iprot),
2022 & maxdimnat*natlike(iprot)*icount,
2023 & MPI_DOUBLE_PRECISION,
2024 & Next,577,WHAM_COMM,IERROR)
2025 call MPI_Recv(eini(jj_old,iprot),icount,
2026 & MPI_DOUBLE_PRECISION,Previous,578,WHAM_COMM,STATUS,IERROR)
2027 call MPI_Send(eini(jj_old,iprot),icount,
2028 & MPI_DOUBLE_PRECISION,Next,578,WHAM_COMM,IERROR)
2029 call MPI_Recv(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2030 & Previous,579,WHAM_COMM,STATUS,IERROR)
2031 call MPI_Send(entfac(jj_old,iprot),icount,MPI_DOUBLE_PRECISION,
2032 & Next,579,WHAM_COMM,IERROR)
2033 call MPI_Recv(c_zs(1,1,1,iprot),3*icount*maxres2,
2034 & MPI_REAL,Previous,580,WHAM_COMM,STATUS,IERROR)
2035 call MPI_Send(c_zs(1,1,1,iprot),3*icount*maxres2,
2036 & MPI_REAL,Next,580,WHAM_COMM,IERROR)
2037 if (.not. init_ene) then
2038 call MPI_Recv(eneps(1,1,1,iprot),2*icount*nntyp,
2039 & MPI_DOUBLE_PRECISION,Previous,582,WHAM_COMM,STATUS,IERROR)
2040 call MPI_Send(eneps(1,1,1,iprot),2*icount*nntyp,
2041 & MPI_DOUBLE_PRECISION,Next,582,WHAM_COMM,IERROR)
2045 list_conf(i,iprot)=i
2047 call dawrite_ccoords(iprot,jj_old,jj,ientout)
2048 call dawrite_ene(iprot,jj_old,jj,istat)
2051 write (iout,*) "Protein",iprot
2052 write (iout,*) "Processor",me," received",icount," conformations"
2054 write (iout,'(8f10.4)') ((c_zs(l,k,i,iprot),l=1,3,k=1,nres)
2055 write (iout,'(8f10.4)')((c_zs(l,k,i+nres,iprot),l=1,3,k=nnt,nct)
2056 write (iout,'(16i5)') nss_zs(i,iprot),(ihpb_zs(k,i,iprot),
2057 & jhpb_zs(k,i,iprot),k=1,nss_zs(i,iprot))
2058 write (iout,'(5e15.5)') (enetb(i,j,iprot),j=1,n_ene)
2059 write (iout,'(2e15.5)') (eneps(1,j,i,iprot),
2060 & eneps(2,j,i,iprot),j=1,nntyp)
2061 write (iout,'(e15.5,16i5)') entfac(i),iscore(i,0,iprot),
2069 c------------------------------------------------------------------------------
2070 subroutine read_thermal
2072 include "DIMENSIONS"
2073 include "DIMENSIONS.ZSCOPT"
2076 include "COMMON.MPI"
2078 include "COMMON.CHAIN"
2079 include "COMMON.SBRIDGE"
2080 include "COMMON.INTERACT"
2081 include "COMMON.IOUNITS"
2082 include "COMMON.CLASSES"
2083 include "COMMON.VAR"
2084 include "COMMON.THERMAL"
2085 character*800 controlcard
2086 call card_concat(controlcard,.true.)
2087 call readi(controlcard,"NGRIDT",NGridT,200)
2088 call reada(controlcard,"DELTAT",deltaT,5.0d0)
2089 call reada(controlcard,"T0",GridT0,2.0d2)
2090 if (me.eq.Master) then
2091 write (iout,*) "Parameters of thermal curves"
2092 write (iout,*) "NGRIDT",NGridT," DELTAT",deltaT," T0",GridT0
2096 c------------------------------------------------------------------------------
2097 subroutine daread_ene(iprot,istart_conf,iend_conf)
2099 include "DIMENSIONS"
2100 include "DIMENSIONS.ZSCOPT"
2103 include "COMMON.MPI"
2105 include "COMMON.CHAIN"
2106 include "COMMON.CLASSES"
2107 include "COMMON.ENERGIES"
2108 include "COMMON.IOUNITS"
2109 include "COMMON.PROTFILES"
2110 include "COMMON.ALLPROT"
2111 include "COMMON.WEIGHTDER"
2112 include "COMMON.COMPAR"
2113 include "COMMON.VMCPAR"
2114 integer iprot,istart_conf,iend_conf
2115 integer i,ii,iii,j,k
2117 write (iout,*) "Calling DAREAD_ENE"
2119 c write (iout,*) "Reading: n_ene",n_ene
2121 c write (iout,*) "DAREAD_ENE",istart_conf,iend_conf
2123 c Read conformations off a DA scratchfile.
2125 do ii=istart_conf,iend_conf
2126 iii=list_conf(ii,iprot)
2127 i = ii - istart_conf + 1
2128 read(ientin,rec=iii) (enetb(i,j,iprot),j=1,n_ene),
2129 & (enetb_orig(i,j,iprot),j=1,n_ene),
2130 & (enetb_oorig(i,j,iprot),j=1,n_ene),
2131 & eini(ii,iprot),entfac(ii,iprot),
2132 & (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp),
2133 & ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2135 write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot),
2137 write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene)
2138 write (iout,'(18(1pe12.4))')
2139 & ((nu(k,j,i,iprot)k=1,maxdimnat),j=1,natlike(iprot))
2146 c------------------------------------------------------------------------------
2147 subroutine dawrite_ene(iprot,istart_conf,iend_conf,unit_out)
2149 include "DIMENSIONS"
2150 include "DIMENSIONS.ZSCOPT"
2153 include "COMMON.MPI"
2155 include "COMMON.CHAIN"
2156 include "COMMON.CLASSES"
2157 include "COMMON.ENERGIES"
2158 include "COMMON.IOUNITS"
2159 include "COMMON.PROTFILES"
2160 include "COMMON.ALLPROT"
2161 include "COMMON.WEIGHTDER"
2162 include "COMMON.VMCPAR"
2163 include "COMMON.COMPAR"
2164 integer iprot,istart_conf,iend_conf,unit_out
2165 integer i,ii,iii,j,k
2166 c write (iout,*) "Writing: n_ene",n_ene
2168 c write (iout,*) "DAWRITE_ENE",istart_conf,iend_conf
2170 c Write conformations to a DA scratchfile.
2172 do ii=istart_conf,iend_conf
2173 iii=list_conf(ii,iprot)
2174 i = ii - istart_conf + 1
2175 write(unit_out,rec=iii) (enetb(i,j,iprot),j=1,n_ene),
2176 & (enetb_orig(i,j,iprot),j=1,n_ene),
2177 & (enetb_oorig(i,j,iprot),j=1,n_ene),
2178 & eini(ii,iprot),entfac(ii,iprot),
2179 & (eneps(1,j,i,iprot),eneps(2,j,i,iprot),j=1,nntyp),
2180 & ((nu(k,j,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2182 write (iout,'(3i5,3e15.4,i5,i10)') ii,iii,i,eini(ii,iprot),
2184 write (iout,'(20(1pe12.4)') (enetb(i,j,iprot),j=1,n_ene)
2185 write (iout,'(18(1pe12.4))')
2186 & ((nu(kj,i,iprot),k=1,maxdimnat),j=1,natlike(iprot))
2193 c------------------------------------------------------------------------------
2194 subroutine daread_ccoords(iprot,istart_conf,iend_conf)
2196 include "DIMENSIONS"
2197 include "DIMENSIONS.ZSCOPT"
2200 include "COMMON.MPI"
2202 include "COMMON.CHAIN"
2203 include "COMMON.CLASSES"
2204 include "COMMON.ENERGIES"
2205 include "COMMON.IOUNITS"
2206 include "COMMON.PROTFILES"
2207 include "COMMON.ALLPROT"
2208 include "COMMON.INTERACT"
2209 include "COMMON.VAR"
2210 include "COMMON.SBRIDGE"
2211 include "COMMON.GEO"
2212 include "COMMON.COMPAR"
2213 include "COMMON.VMCPAR"
2214 include "COMMON.WEIGHTDER"
2215 integer iprot,istart_conf,iend_conf
2216 integer i,j,k,ij,ii,iii
2218 real*4 csingle(3,maxres2)
2219 character*16 form,acc
2222 c Read conformations off a DA scratchfile.
2225 write (iout,*) "DAREAD_COORDS"
2226 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2227 write (iout,*) "lenrec",lenrec(iprot)
2228 inquire(unit=ientin,name=nam,recl=len,form=form,access=acc)
2229 write (iout,*) "len=",len," form=",form," acc=",acc
2230 write (iout,*) "nam=",nam
2233 do ii=istart_conf,iend_conf
2234 iii=list_conf(ii,iprot)
2235 ij = ii - istart_conf + 1
2237 write (iout,*) "Reading binary file, record",iii," ii",ii
2240 read(icbase,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2241 & ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2242 & nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot),
2244 & ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot))
2251 write (iout,*) "iprot",iprot," ii",ii
2252 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2253 write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2254 write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot)
2255 write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2258 call store_ccoords(iprot,ii-istart_conf+1)
2262 c------------------------------------------------------------------------------
2263 subroutine dawrite_ccoords(iprot,istart_conf,iend_conf,unit_out)
2265 include "DIMENSIONS"
2266 include "DIMENSIONS.ZSCOPT"
2269 include "COMMON.MPI"
2271 include "COMMON.CHAIN"
2272 include "COMMON.INTERACT"
2273 include "COMMON.CLASSES"
2274 include "COMMON.ENERGIES"
2275 include "COMMON.IOUNITS"
2276 include "COMMON.PROTFILES"
2277 include "COMMON.ALLPROT"
2278 include "COMMON.VAR"
2279 include "COMMON.SBRIDGE"
2280 include "COMMON.GEO"
2281 include "COMMON.COMPAR"
2282 include "COMMON.WEIGHTDER"
2283 include "COMMON.VMCPAR"
2284 real*4 csingle(3,maxres2)
2285 integer iprot,istart_conf,iend_conf
2286 integer i,j,k,ii,ij,iii,unit_out
2288 character*16 form,acc
2291 c Write conformations to a DA scratchfile.
2294 write (iout,*) "DAWRITE_COORDS"
2295 write (iout,*) "istart_conf",istart_conf," iend_conf",iend_conf
2296 write (iout,*) "lenrec",lenrec(iprot)
2297 inquire(unit=unit_out,name=nam,recl=len,form=form,access=acc)
2298 write (iout,*) "len=",len," form=",form," acc=",acc
2299 write (iout,*) "nam=",nam
2302 do ii=istart_conf,iend_conf
2303 iii=list_conf(ii,iprot)
2304 ij = ii - istart_conf + 1
2305 call restore_ccoords(iprot,ii-istart_conf+1)
2307 write (iout,*) "Writing binary file, record",iii," ii",ii
2315 write(unit_out,rec=iii) ((csingle(j,i),j=1,3),i=1,nres),
2316 & ((csingle(j,i),j=1,3),i=nnt+nres,nct+nres),
2317 & nss,(ihpb(i),jhpb(i),i=1,nss),eini(ii,iprot),
2319 & ((nu(k,i,ij,iprot),k=1,maxdimnat),i=1,natlike(iprot))
2321 write (iout,*) "kbatch",kbatch(ii,iprot)," iscore",
2322 & iscore(ii,0,iprot)
2323 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2324 write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2325 write (iout,'(2e15.5)') eini(ii,iprot),entfac(ii,iprot)
2326 write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2332 c------------------------------------------------------------------------------
2333 subroutine store_ccoords(iprot,ii)
2335 include "DIMENSIONS"
2336 include "DIMENSIONS.ZSCOPT"
2337 include "COMMON.VAR"
2338 include "COMMON.CHAIN"
2339 include "COMMON.ALLPROT"
2340 include "COMMON.IOUNITS"
2341 include "COMMON.GEO"
2342 include "COMMON.SBRIDGE"
2343 double precision thetnorm
2344 integer iprot,i,j,ii
2345 do i=1,nres_zs(iprot)
2347 c_zs(j,i,ii,iprot)=c(j,i)
2350 do i=nnt_zs(iprot),nct_zs(iprot)
2352 c_zs(j,i+nres,ii,iprot)=c(j,i+nres)
2355 c 5/7/02 AL: store sbridge info
2356 nss_zs(ii,iprot)=nss
2358 ihpb_zs(i,ii,iprot)=ihpb(i)
2359 jhpb_zs(i,ii,iprot)=jhpb(i)
2363 c------------------------------------------------------------------------------
2364 subroutine restore_ccoords(iprot,ii)
2366 include "DIMENSIONS"
2367 include "DIMENSIONS.ZSCOPT"
2368 include "COMMON.INTERACT"
2369 include "COMMON.VAR"
2370 include "COMMON.ALLPROT"
2371 include "COMMON.SBRIDGE"
2372 include "COMMON.CHAIN"
2373 include "COMMON.IOUNITS"
2374 integer iprot,i,j,ii
2375 do i=1,nres_zs(iprot)
2377 c(j,i)=c_zs(j,i,ii,iprot)
2380 do i=nnt_zs(iprot),nct_zs(iprot)
2382 c(j,i+nres)=c_zs(j,i+nres,ii,iprot)
2385 c 5/7/02 AL: restore sbridge info
2386 nss=nss_zs(ii,iprot)
2388 ihpb(i)=ihpb_zs(i,ii,iprot)+nres
2389 jhpb(i)=jhpb_zs(i,ii,iprot)+nres
2394 write (iout,*) "restore_ccoords",ii
2395 write (iout,'(8f10.5)') ((c(j,i),j=1,3),i=1,nres)
2396 write (iout,'(8f10.4)') ((c(j,i),j=1,3),i=nnt+nres,nct+nres)
2397 write (iout,'(16i5)') nss,(ihpb(i),jhpb(i),i=1,nss)
2402 c------------------------------------------------------------------------------
2403 subroutine card_concat(card,to_upper)
2405 include 'DIMENSIONS.ZSCOPT'
2406 include "COMMON.IOUNITS"
2408 character*80 karta,ucase
2412 read (inp,'(a)') karta
2413 if (to_upper) karta=ucase(karta)
2415 do while (karta(80:80).eq.'&')
2416 card=card(:ilen(card)+1)//karta(:79)
2417 read (inp,'(a)') karta
2418 if (to_upper) karta=ucase(karta)
2420 card=card(:ilen(card)+1)//karta
2423 c------------------------------------------------------------------------------
2424 subroutine readi(rekord,lancuch,wartosc,default)
2426 character*(*) rekord,lancuch
2427 integer wartosc,default
2430 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2431 if (iread.eq.0) then
2435 iread=iread+ilen(lancuch)+1
2436 read (rekord(iread:),*) wartosc
2439 c----------------------------------------------------------------------------
2440 subroutine reada(rekord,lancuch,wartosc,default)
2442 character*(*) rekord,lancuch
2444 double precision wartosc,default
2447 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2448 if (iread.eq.0) then
2452 iread=iread+ilen(lancuch)+1
2453 read (rekord(iread:),*) wartosc
2456 c----------------------------------------------------------------------------
2457 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2460 integer tablica(dim),default
2461 character*(*) rekord,lancuch
2468 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2469 if (iread.eq.0) return
2470 iread=iread+ilen(lancuch)+1
2471 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2474 c----------------------------------------------------------------------------
2475 subroutine multreada(rekord,lancuch,tablica,dim,default)
2478 double precision tablica(dim),default
2479 character*(*) rekord,lancuch
2486 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2487 if (iread.eq.0) return
2488 iread=iread+ilen(lancuch)+1
2489 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2492 c----------------------------------------------------------------------------
2493 subroutine reads(rekord,lancuch,wartosc,default)
2495 character*(*) rekord,lancuch,wartosc,default
2497 integer ilen,lenlan,lenrec,iread,ireade
2501 lenlan=ilen(lancuch)
2503 iread=index(rekord,lancuch(:lenlan)//"=")
2504 c print *,"rekord",rekord," lancuch",lancuch
2505 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2506 if (iread.eq.0) then
2510 iread=iread+lenlan+1
2511 c print *,"iread",iread
2512 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2513 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2515 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2517 c print *,"iread",iread
2518 if (iread.gt.lenrec) then
2523 c print *,"ireade",ireade
2524 do while (ireade.lt.lenrec .and.
2525 & .not.iblnk(rekord(ireade:ireade)))
2528 wartosc=rekord(iread:ireade)
2531 c----------------------------------------------------------------------------
2532 subroutine multreads(rekord,lancuch,tablica,dim,default)
2535 character*(*) rekord,lancuch,tablica(dim),default
2537 integer ilen,lenlan,lenrec,iread,ireade
2544 lenlan=ilen(lancuch)
2546 iread=index(rekord,lancuch(:lenlan)//"=")
2547 c print *,"rekord",rekord," lancuch",lancuch
2548 c print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
2549 if (iread.eq.0) return
2550 iread=iread+lenlan+1
2552 c print *,"iread",iread
2553 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2554 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
2556 c print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
2558 c print *,"iread",iread
2559 if (iread.gt.lenrec) return
2561 c print *,"ireade",ireade
2562 do while (ireade.lt.lenrec .and.
2563 & .not.iblnk(rekord(ireade:ireade)))
2566 tablica(i)=rekord(iread:ireade)
2570 c----------------------------------------------------------------------------
2571 subroutine split_string(rekord,tablica,dim,nsub)
2573 integer dim,nsub,i,ii,ll,kk
2574 character*(*) tablica(dim)
2575 character*(*) rekord
2585 C Find the start of term name
2587 do while (ii.le.ll .and. rekord(ii:ii).eq." ")
2590 C Parse the name into TABLICA(i) until blank found
2591 do while (ii.le.ll .and. rekord(ii:ii).ne." ")
2593 tablica(i)(kk:kk)=rekord(ii:ii)
2596 if (kk.gt.0) nsub=nsub+1
2597 if (ii.gt.ll) return
2601 c-------------------------------------------------------------------------------
2602 integer function iroof(n,m)
2604 if (ii*m .lt. n) ii=ii+1