2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 write (iout,'(/a,i5,a/4a5,2a8,3a10,a5)')
46 & "The following",nhpb-nss,
47 & " distance restraints have been imposed:",
48 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
51 write (iout,'(4i5,2f8.2,3f10.5,i5)')i-nss,ihpb(i),jhpb(i),
52 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
57 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
58 & "The following",npeak,
59 & " NMR peak restraints have been imposed:",
60 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
61 & " score"," type"," ipeak"
63 do j=ipeak(1,i),ipeak(2,i)
64 write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
65 & jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
66 & forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
69 write (iout,*) "The ipeak array"
71 write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
77 c print *,"Processor",myrank," leaves READRTNS"
80 C-------------------------------------------------------------------------------
81 subroutine read_control
85 implicit real*8 (a-h,o-z)
89 logical OKRandom, prng_restart
92 include 'COMMON.IOUNITS'
93 include 'COMMON.TIME1'
94 include 'COMMON.THREAD'
95 include 'COMMON.SBRIDGE'
96 include 'COMMON.CONTROL'
99 include 'COMMON.HEADER'
101 include 'COMMON.CHAIN'
102 include 'COMMON.MUCA'
104 include 'COMMON.FFIELD'
105 include 'COMMON.INTERACT'
106 include 'COMMON.SETUP'
107 include 'COMMON.SPLITELE'
108 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
109 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
111 character*320 controlcard
117 read (INP,'(a)') titel
118 call card_concat(controlcard)
119 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
120 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
121 call reada(controlcard,'SEED',seed,0.0D0)
122 call random_init(seed)
123 C Set up the time limit (caution! The time must be input in minutes!)
124 read_cart=index(controlcard,'READ_CART').gt.0
125 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
126 write (iout,*) "constr_dist",constr_dist
127 call readi(controlcard,'NSAXS',nsaxs,0)
128 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
129 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
130 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
131 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
132 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
133 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
134 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
135 call readi(controlcard,'SYM',symetr,1)
136 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
137 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
138 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
139 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
140 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
141 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
142 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
143 call reada(controlcard,'DRMS',drms,0.1D0)
144 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
145 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
146 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
147 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
148 write (iout,'(a,f10.1)')'DRMS = ',drms
149 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
150 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
152 call readi(controlcard,'NZ_START',nz_start,0)
153 call readi(controlcard,'NZ_END',nz_end,0)
154 call readi(controlcard,'IZ_SC',iz_sc,0)
156 safety = 60.0d0*safety
159 call reada(controlcard,"T_BATH",t_bath,300.0d0)
160 minim=(index(controlcard,'MINIMIZE').gt.0)
161 dccart=(index(controlcard,'CART').gt.0)
162 overlapsc=(index(controlcard,'OVERLAP').gt.0)
163 overlapsc=.not.overlapsc
164 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
165 searchsc=.not.searchsc
166 sideadd=(index(controlcard,'SIDEADD').gt.0)
167 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
168 outpdb=(index(controlcard,'PDBOUT').gt.0)
169 outmol2=(index(controlcard,'MOL2OUT').gt.0)
170 pdbref=(index(controlcard,'PDBREF').gt.0)
171 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
172 indpdb=index(controlcard,'PDBSTART')
173 extconf=(index(controlcard,'EXTCONF').gt.0)
174 AFMlog=(index(controlcard,'AFM'))
175 selfguide=(index(controlcard,'SELFGUIDE'))
176 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
177 call readi(controlcard,'IPRINT',iprint,0)
178 call readi(controlcard,'MAXGEN',maxgen,10000)
179 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
180 call readi(controlcard,"KDIAG",kdiag,0)
181 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
182 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
183 & write (iout,*) "RESCALE_MODE",rescale_mode
184 split_ene=index(controlcard,'SPLIT_ENE').gt.0
185 if (index(controlcard,'REGULAR').gt.0.0D0) then
186 call reada(controlcard,'WEIDIS',weidis,0.1D0)
190 if (index(controlcard,'CHECKGRAD').gt.0) then
192 if (index(controlcard,'CART').gt.0) then
194 elseif (index(controlcard,'CARINT').gt.0) then
199 elseif (index(controlcard,'THREAD').gt.0) then
201 call readi(controlcard,'THREAD',nthread,0)
202 if (nthread.gt.0) then
203 call reada(controlcard,'WEIDIS',weidis,0.1D0)
206 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
207 stop 'Error termination in Read_Control.'
209 else if (index(controlcard,'MCMA').gt.0) then
211 else if (index(controlcard,'MCEE').gt.0) then
213 else if (index(controlcard,'MULTCONF').gt.0) then
215 else if (index(controlcard,'MAP').gt.0) then
217 call readi(controlcard,'MAP',nmap,0)
218 else if (index(controlcard,'CSA').gt.0) then
220 crc else if (index(controlcard,'ZSCORE').gt.0) then
222 crc ZSCORE is rm from UNRES, modecalc=9 is available
225 cfcm else if (index(controlcard,'MCMF').gt.0) then
227 else if (index(controlcard,'SOFTREG').gt.0) then
229 else if (index(controlcard,'CHECK_BOND').gt.0) then
231 else if (index(controlcard,'TEST').gt.0) then
233 else if (index(controlcard,'MD').gt.0) then
235 else if (index(controlcard,'RE ').gt.0) then
239 lmuca=index(controlcard,'MUCA').gt.0
240 call readi(controlcard,'MUCADYN',mucadyn,0)
241 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
242 if (lmuca .and. (me.eq.king .or. .not.out1file ))
244 write (iout,*) 'MUCADYN=',mucadyn
245 write (iout,*) 'MUCASMOOTH=',muca_smooth
248 iscode=index(controlcard,'ONE_LETTER')
249 indphi=index(controlcard,'PHI')
250 indback=index(controlcard,'BACK')
251 iranconf=index(controlcard,'RAND_CONF')
252 i2ndstr=index(controlcard,'USE_SEC_PRED')
253 gradout=index(controlcard,'GRADOUT').gt.0
254 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
255 C DISTCHAINMAX become obsolete for periodic boundry condition
256 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
257 C Reading the dimensions of box in x,y,z coordinates
258 call reada(controlcard,'BOXX',boxxsize,100.0d0)
259 call reada(controlcard,'BOXY',boxysize,100.0d0)
260 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
261 c Cutoff range for interactions
262 call reada(controlcard,"R_CUT",r_cut,15.0d0)
263 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
264 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
265 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
266 if (lipthick.gt.0.0d0) then
267 bordliptop=(boxzsize+lipthick)/2.0
268 bordlipbot=bordliptop-lipthick
270 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
271 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
272 buflipbot=bordlipbot+lipbufthick
273 bufliptop=bordliptop-lipbufthick
274 if ((lipbufthick*2.0d0).gt.lipthick)
275 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
277 c write(iout,*) "bordliptop=",bordliptop
278 c write(iout,*) "bordlipbot=",bordlipbot
279 c write(iout,*) "bufliptop=",bufliptop
280 c write(iout,*) "buflipbot=",buflipbot
283 if (me.eq.king .or. .not.out1file )
284 & write (iout,*) "DISTCHAINMAX",distchainmax
286 if(me.eq.king.or..not.out1file)
287 & write (iout,'(2a)') diagmeth(kdiag),
288 & ' routine used to diagonalize matrices.'
291 c--------------------------------------------------------------------------
292 subroutine read_REMDpar
296 implicit real*8 (a-h,o-z)
298 include 'COMMON.IOUNITS'
299 include 'COMMON.TIME1'
302 include 'COMMON.LANGEVIN'
304 include 'COMMON.LANGEVIN.lang0'
306 include 'COMMON.INTERACT'
307 include 'COMMON.NAMES'
309 include 'COMMON.REMD'
310 include 'COMMON.CONTROL'
311 include 'COMMON.SETUP'
313 character*320 controlcard
314 character*3200 controlcard1
315 integer iremd_m_total
317 if(me.eq.king.or..not.out1file)
318 & write (iout,*) "REMD setup"
320 call card_concat(controlcard)
321 call readi(controlcard,"NREP",nrep,3)
322 call readi(controlcard,"NSTEX",nstex,1000)
323 call reada(controlcard,"RETMIN",retmin,10.0d0)
324 call reada(controlcard,"RETMAX",retmax,1000.0d0)
325 mremdsync=(index(controlcard,'SYNC').gt.0)
326 call readi(controlcard,"NSYN",i_sync_step,100)
327 restart1file=(index(controlcard,'REST1FILE').gt.0)
328 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
329 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
330 if(max_cache_traj_use.gt.max_cache_traj)
331 & max_cache_traj_use=max_cache_traj
332 if(me.eq.king.or..not.out1file) then
333 cd if (traj1file) then
334 crc caching is in testing - NTWX is not ignored
335 cd write (iout,*) "NTWX value is ignored"
336 cd write (iout,*) " trajectory is stored to one file by master"
337 cd write (iout,*) " before exchange at NSTEX intervals"
339 write (iout,*) "NREP= ",nrep
340 write (iout,*) "NSTEX= ",nstex
341 write (iout,*) "SYNC= ",mremdsync
342 write (iout,*) "NSYN= ",i_sync_step
343 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
346 if (index(controlcard,'TLIST').gt.0) then
348 call card_concat(controlcard1)
349 read(controlcard1,*) (remd_t(i),i=1,nrep)
350 if(me.eq.king.or..not.out1file)
351 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
354 if (index(controlcard,'MLIST').gt.0) then
356 call card_concat(controlcard1)
357 read(controlcard1,*) (remd_m(i),i=1,nrep)
358 if(me.eq.king.or..not.out1file) then
359 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
362 iremd_m_total=iremd_m_total+remd_m(i)
364 write (iout,*) 'Total number of replicas ',iremd_m_total
367 if(me.eq.king.or..not.out1file)
368 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
371 c--------------------------------------------------------------------------
372 subroutine read_MDpar
376 implicit real*8 (a-h,o-z)
378 include 'COMMON.IOUNITS'
379 include 'COMMON.TIME1'
382 include 'COMMON.LANGEVIN'
384 include 'COMMON.LANGEVIN.lang0'
386 include 'COMMON.INTERACT'
387 include 'COMMON.NAMES'
389 include 'COMMON.SETUP'
390 include 'COMMON.CONTROL'
391 include 'COMMON.SPLITELE'
393 character*320 controlcard
395 call card_concat(controlcard)
396 call readi(controlcard,"NSTEP",n_timestep,1000000)
397 call readi(controlcard,"NTWE",ntwe,100)
398 call readi(controlcard,"NTWX",ntwx,1000)
399 call reada(controlcard,"DT",d_time,1.0d-1)
400 call reada(controlcard,"DVMAX",dvmax,2.0d1)
401 call reada(controlcard,"DAMAX",damax,1.0d1)
402 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
403 call readi(controlcard,"LANG",lang,0)
404 RESPA = index(controlcard,"RESPA") .gt. 0
405 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
406 ntime_split0=ntime_split
407 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
408 ntime_split0=ntime_split
409 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
410 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
411 rest = index(controlcard,"REST").gt.0
412 tbf = index(controlcard,"TBF").gt.0
413 usampl = index(controlcard,"USAMPL").gt.0
415 mdpdb = index(controlcard,"MDPDB").gt.0
416 call reada(controlcard,"T_BATH",t_bath,300.0d0)
417 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
418 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
419 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
420 if (count_reset_moment.eq.0) count_reset_moment=1000000000
421 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
422 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
423 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
424 if (count_reset_vel.eq.0) count_reset_vel=1000000000
425 large = index(controlcard,"LARGE").gt.0
426 print_compon = index(controlcard,"PRINT_COMPON").gt.0
427 rattle = index(controlcard,"RATTLE").gt.0
428 preminim = index(controlcard,"PREMINIM").gt.0
430 dccart=(index(controlcard,'CART').gt.0)
433 c if performing umbrella sampling, fragments constrained are read from the fragment file
439 if(me.eq.king.or..not.out1file) then
441 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
443 write (iout,'(a)') "The units are:"
444 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
445 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
446 & " acceleration: angstrom/(48.9 fs)**2"
447 write (iout,'(a)') "energy: kcal/mol, temperature: K"
449 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
450 write (iout,'(a60,f10.5,a)')
451 & "Initial time step of numerical integration:",d_time,
453 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
455 write (iout,'(2a,i4,a)')
456 & "A-MTS algorithm used; initial time step for fast-varying",
457 & " short-range forces split into",ntime_split," steps."
458 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
459 & r_cut," lambda",rlamb
461 write (iout,'(2a,f10.5)')
462 & "Maximum acceleration threshold to reduce the time step",
463 & "/increase split number:",damax
464 write (iout,'(2a,f10.5)')
465 & "Maximum predicted energy drift to reduce the timestep",
466 & "/increase split number:",edriftmax
467 write (iout,'(a60,f10.5)')
468 & "Maximum velocity threshold to reduce velocities:",dvmax
469 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
470 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
471 if (rattle) write (iout,'(a60)')
472 & "Rattle algorithm used to constrain the virtual bonds"
473 if (preminim .or. iranconf.gt.0) then
475 & "Initial structure will be energy-minimized"
480 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
481 call reada(controlcard,"RWAT",rwat,1.4d0)
482 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
483 surfarea=index(controlcard,"SURFAREA").gt.0
484 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
485 if(me.eq.king.or..not.out1file)then
486 write (iout,'(/a,$)') "Langevin dynamics calculation"
489 & " with direct integration of Langevin equations"
490 else if (lang.eq.2) then
491 write (iout,'(a/)') " with TINKER stochasic MD integrator"
492 else if (lang.eq.3) then
493 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
494 else if (lang.eq.4) then
495 write (iout,'(a/)') " in overdamped mode"
497 write (iout,'(//a,i5)')
498 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
501 write (iout,'(a60,f10.5)') "Temperature:",t_bath
502 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
503 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
504 write (iout,'(a60,f10.5)')
505 & "Scaling factor of the friction forces:",scal_fric
506 if (surfarea) write (iout,'(2a,i10,a)')
507 & "Friction coefficients will be scaled by solvent-accessible",
508 & " surface area every",reset_fricmat," steps."
510 c Calculate friction coefficients and bounds of stochastic forces
511 eta=6*pi*cPoise*etawat
512 if(me.eq.king.or..not.out1file)
513 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
515 gamp=scal_fric*(pstok+rwat)*eta
516 stdfp=dsqrt(2*Rb*t_bath/d_time)
518 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
519 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
521 if(me.eq.king.or..not.out1file)then
522 write (iout,'(/2a/)')
523 & "Radii of site types and friction coefficients and std's of",
524 & " stochastic forces of fully exposed sites"
525 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
527 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
528 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
532 if(me.eq.king.or..not.out1file)then
533 write (iout,'(a)') "Berendsen bath calculation"
534 write (iout,'(a60,f10.5)') "Temperature:",t_bath
535 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
537 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
538 & count_reset_moment," steps"
540 & write (iout,'(a,i10,a)')
541 & "Velocities will be reset at random every",count_reset_vel,
545 if(me.eq.king.or..not.out1file)
546 & write (iout,'(a31)') "Microcanonical mode calculation"
548 if(me.eq.king.or..not.out1file)then
549 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
551 write(iout,*) "MD running with constraints."
552 write(iout,*) "Equilibration time ", eq_time, " mtus."
553 write(iout,*) "Constraining ", nfrag," fragments."
554 write(iout,*) "Length of each fragment, weight and q0:"
556 write (iout,*) "Set of restraints #",iset
558 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
559 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
561 write(iout,*) "constraints between ", npair, "fragments."
562 write(iout,*) "constraint pairs, weights and q0:"
564 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
565 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
567 write(iout,*) "angle constraints within ", nfrag_back,
568 & "backbone fragments."
569 write(iout,*) "fragment, weights:"
571 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
572 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
573 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
576 iset=mod(kolor,nset)+1
579 if(me.eq.king.or..not.out1file)
580 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
583 c------------------------------------------------------------------------------
586 C Read molecular data.
588 implicit real*8 (a-h,o-z)
594 include 'COMMON.IOUNITS'
597 include 'COMMON.INTERACT'
598 include 'COMMON.LOCAL'
599 include 'COMMON.NAMES'
600 include 'COMMON.CHAIN'
601 include 'COMMON.FFIELD'
602 include 'COMMON.SBRIDGE'
603 include 'COMMON.HEADER'
604 include 'COMMON.CONTROL'
605 include 'COMMON.DBASE'
606 include 'COMMON.THREAD'
607 include 'COMMON.CONTACTS'
608 include 'COMMON.TORCNSTR'
609 include 'COMMON.TIME1'
610 include 'COMMON.BOUNDS'
612 include 'COMMON.SETUP'
613 character*4 sequence(maxres)
615 double precision x(maxvar)
616 character*256 pdbfile
617 character*320 weightcard
618 character*80 weightcard_t,ucase
619 dimension itype_pdb(maxres)
620 common /pizda/ itype_pdb
621 logical seq_comp,fail
622 double precision energia(0:n_ene)
628 C Read weights of the subsequent energy terms.
629 call card_concat(weightcard)
630 call reada(weightcard,'WLONG',wlong,1.0D0)
631 call reada(weightcard,'WSC',wsc,wlong)
632 call reada(weightcard,'WSCP',wscp,wlong)
633 call reada(weightcard,'WELEC',welec,1.0D0)
634 call reada(weightcard,'WVDWPP',wvdwpp,welec)
635 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
636 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
637 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
638 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
639 call reada(weightcard,'WTURN3',wturn3,1.0D0)
640 call reada(weightcard,'WTURN4',wturn4,1.0D0)
641 call reada(weightcard,'WTURN6',wturn6,1.0D0)
642 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
643 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
644 call reada(weightcard,'WBOND',wbond,1.0D0)
645 call reada(weightcard,'WTOR',wtor,1.0D0)
646 call reada(weightcard,'WTORD',wtor_d,1.0D0)
647 call reada(weightcard,'WANG',wang,1.0D0)
648 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
649 call reada(weightcard,'SCAL14',scal14,0.4D0)
650 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
651 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
652 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
653 call reada(weightcard,'TEMP0',temp0,300.0d0)
654 call reada(weightcard,'WLT',wliptran,0.0D0)
655 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
656 if (index(weightcard,'SOFT').gt.0) ipot=6
657 C 12/1/95 Added weight for the multi-body term WCORR
658 call reada(weightcard,'WCORRH',wcorr,1.0D0)
659 if (wcorr4.gt.0.0d0) wcorr=wcorr4
680 if(me.eq.king.or..not.out1file)
681 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
682 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
684 10 format (/'Energy-term weights (unscaled):'//
685 & 'WSCC= ',f10.6,' (SC-SC)'/
686 & 'WSCP= ',f10.6,' (SC-p)'/
687 & 'WELEC= ',f10.6,' (p-p electr)'/
688 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
689 & 'WBOND= ',f10.6,' (stretching)'/
690 & 'WANG= ',f10.6,' (bending)'/
691 & 'WSCLOC= ',f10.6,' (SC local)'/
692 & 'WTOR= ',f10.6,' (torsional)'/
693 & 'WTORD= ',f10.6,' (double torsional)'/
694 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
695 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
696 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
697 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
698 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
699 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
700 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
701 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
702 & 'WTURN6= ',f10.6,' (turns, 6th order)')
703 if(me.eq.king.or..not.out1file)then
704 if (wcorr4.gt.0.0d0) then
705 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
706 & 'between contact pairs of peptide groups'
707 write (iout,'(2(a,f5.3/))')
708 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
709 & 'Range of quenching the correlation terms:',2*delt_corr
710 else if (wcorr.gt.0.0d0) then
711 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
712 & 'between contact pairs of peptide groups'
714 write (iout,'(a,f8.3)')
715 & 'Scaling factor of 1,4 SC-p interactions:',scal14
716 write (iout,'(a,f8.3)')
717 & 'General scaling factor of SC-p interactions:',scalscp
719 r0_corr=cutoff_corr-delt_corr
721 aad(i,1)=scalscp*aad(i,1)
722 aad(i,2)=scalscp*aad(i,2)
723 bad(i,1)=scalscp*bad(i,1)
724 bad(i,2)=scalscp*bad(i,2)
726 call rescale_weights(t_bath)
727 if(me.eq.king.or..not.out1file)
728 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
729 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
731 22 format (/'Energy-term weights (scaled):'//
732 & 'WSCC= ',f10.6,' (SC-SC)'/
733 & 'WSCP= ',f10.6,' (SC-p)'/
734 & 'WELEC= ',f10.6,' (p-p electr)'/
735 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
736 & 'WBOND= ',f10.6,' (stretching)'/
737 & 'WANG= ',f10.6,' (bending)'/
738 & 'WSCLOC= ',f10.6,' (SC local)'/
739 & 'WTOR= ',f10.6,' (torsional)'/
740 & 'WTORD= ',f10.6,' (double torsional)'/
741 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
742 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
743 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
744 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
745 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
746 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
747 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
748 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
749 & 'WTURN6= ',f10.6,' (turns, 6th order)')
750 if(me.eq.king.or..not.out1file)
751 & write (iout,*) "Reference temperature for weights calculation:",
753 call reada(weightcard,"D0CM",d0cm,3.78d0)
754 call reada(weightcard,"AKCM",akcm,15.1d0)
755 call reada(weightcard,"AKTH",akth,11.0d0)
756 call reada(weightcard,"AKCT",akct,12.0d0)
757 call reada(weightcard,"V1SS",v1ss,-1.08d0)
758 call reada(weightcard,"V2SS",v2ss,7.61d0)
759 call reada(weightcard,"V3SS",v3ss,13.7d0)
760 call reada(weightcard,"EBR",ebr,-5.50D0)
761 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
763 dyn_ss_mask(i)=.false.
767 dyn_ssbond_ij(i,j)=1.0d300
770 call reada(weightcard,"HT",Ht,0.0D0)
772 ss_depth=ebr/wsc-0.25*eps(1,1)
773 Ht=Ht/wsc-0.25*eps(1,1)
774 akcm=akcm*wstrain/wsc
775 akth=akth*wstrain/wsc
776 akct=akct*wstrain/wsc
777 v1ss=v1ss*wstrain/wsc
778 v2ss=v2ss*wstrain/wsc
779 v3ss=v3ss*wstrain/wsc
781 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
784 if(me.eq.king.or..not.out1file) then
785 write (iout,*) "Parameters of the SS-bond potential:"
786 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
788 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
789 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
790 write (iout,*)" HT",Ht
791 print *,'indpdb=',indpdb,' pdbref=',pdbref
793 if (indpdb.gt.0 .or. pdbref) then
794 read(inp,'(a)') pdbfile
795 if(me.eq.king.or..not.out1file)
796 & write (iout,'(2a)') 'PDB data will be read from file ',
797 & pdbfile(:ilen(pdbfile))
798 open(ipdbin,file=pdbfile,status='old',err=33)
800 33 write (iout,'(a)') 'Error opening PDB file.'
803 c print *,'Begin reading pdb data'
812 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
813 & (crefjlee(j,i+nres),j=1,3)
816 c print *,'Finished reading pdb data'
817 if(me.eq.king.or..not.out1file)
818 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
819 & ' nstart_sup=',nstart_sup
821 itype_pdb(i)=itype(i)
825 nct=nstart_sup+nsup-1
826 call contact(.false.,ncont_ref,icont_ref,co)
829 if(me.eq.king.or..not.out1file)
830 & write(iout,*)'Adding sidechains'
834 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
837 do while (fail.and.nsi.le.maxsi)
838 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
842 c Calculalte only the coordinates of the current sidechain; no need to rebuild
844 call locate_side_chain(i)
845 if(fail) write(iout,*)'Adding sidechain failed for res ',
846 & i,' after ',nsi,' trials'
849 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
853 if (indpdb.eq.0) then
854 C Read sequence if not taken from the pdb file.
856 c print *,'nres=',nres
857 if (iscode.gt.0) then
858 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
860 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
862 C Convert sequence to numeric code
864 itype(i)=rescode(i,sequence(i),iscode)
866 C Assign initial virtual bond lengths
872 vbld(i+nres)=dsc(iabs(itype(i)))
873 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
874 c write (iout,*) "i",i," itype",itype(i),
875 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
879 c print '(20i4)',(itype(i),i=1,nres)
882 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
884 if (itype(i).eq.ntyp1) then
888 else if (iabs(itype(i+1)).ne.20) then
890 else if (iabs(itype(i)).ne.20) then
897 if(me.eq.king.or..not.out1file)then
898 write (iout,*) "ITEL"
900 write (iout,*) i,itype(i),itel(i)
902 print *,'Call Read_Bridge.'
905 C 8/13/98 Set limits to generating the dihedral angles
910 read (inp,*) ndih_constr
911 write (iout,*) "ndish_constr",ndih_constr
912 if (ndih_constr.gt.0) then
914 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
915 if(me.eq.king.or..not.out1file)then
917 & 'There are',ndih_constr,' constraints on phi angles.'
919 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
923 phi0(i)=deg2rad*phi0(i)
924 drange(i)=deg2rad*drange(i)
926 if(me.eq.king.or..not.out1file)
927 & write (iout,*) 'FTORS',ftors
930 phibound(1,ii) = phi0(i)-drange(i)
931 phibound(2,ii) = phi0(i)+drange(i)
938 write (iout,'(a)') 'Boundaries in phi angle sampling:'
940 write (iout,'(a3,i5,2f10.1)')
941 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
947 cd print *,'NNT=',NNT,' NCT=',NCT
948 if (itype(1).eq.ntyp1) nnt=2
949 if (itype(nres).eq.ntyp1) nct=nct-1
951 if(me.eq.king.or..not.out1file)
952 & write (iout,'(a,i3)') 'nsup=',nsup
954 if (nsup.le.(nct-nnt+1)) then
955 do i=0,nct-nnt+1-nsup
956 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
962 & 'Error - sequences to be superposed do not match.'
965 do i=0,nsup-(nct-nnt+1)
966 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
968 nstart_sup=nstart_sup+i
974 & 'Error - sequences to be superposed do not match.'
977 if (nsup.eq.0) nsup=nct-nnt
978 if (nstart_sup.eq.0) nstart_sup=nnt
979 if (nstart_seq.eq.0) nstart_seq=nnt
980 if(me.eq.king.or..not.out1file)
981 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
982 & ' nstart_seq=',nstart_seq
984 c--- Zscore rms -------
985 if (nz_start.eq.0) nz_start=nnt
986 if (nz_end.eq.0 .and. nsup.gt.0) then
988 else if (nz_end.eq.0) then
991 if(me.eq.king.or..not.out1file)then
992 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
993 write (iout,*) 'IZ_SC=',iz_sc
995 c----------------------
998 if (.not.pdbref) then
999 call read_angles(inp,*38)
1001 38 write (iout,'(a)') 'Error reading reference structure.'
1003 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1004 stop 'Error reading reference structure'
1006 39 call chainbuild_extconf
1008 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1015 cref(j,i,kkk)=c(j,i)
1018 call contact(.true.,ncont_ref,icont_ref,co)
1022 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1024 if (constr_dist.gt.0) call read_dist_constr
1025 write (iout,*) "After read_dist_constr nhpb",nhpb
1026 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1027 if(me.eq.king.or..not.out1file)
1028 & write (iout,*) 'Contact order:',co
1030 if(me.eq.king.or..not.out1file)
1031 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1034 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1036 if(me.eq.king.or..not.out1file)
1037 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1038 & icont_ref(1,i),' ',
1039 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1042 write (iout,*) "calling read_saxs_consrtr",nsaxs
1043 if (nsaxs.gt.0) call read_saxs_constr
1045 if (constr_homology.gt.0) then
1046 call read_constr_homology
1047 if (indpdb.gt.0 .or. pdbref) then
1050 c(j,i)=crefjlee(j,i)
1051 cref(j,i,1)=crefjlee(j,i)
1056 write (iout,*) "Array C"
1058 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1059 & (c(j,i+nres),j=1,3)
1061 write (iout,*) "Array Cref"
1063 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
1064 & (cref(j,i+nres,1),j=1,3)
1067 call int_from_cart1(.false.)
1068 call sc_loc_geom(.false.)
1070 thetaref(i)=theta(i)
1075 dc(j,i)=c(j,i+1)-c(j,i)
1076 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1081 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1082 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1090 if (nhpb.gt.0) call hpb_partition
1091 if (peak.gt.0) call NMRpeak_partition
1092 c write (iout,*) "After read_dist_constr nhpb",nhpb
1094 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1095 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1096 & modecalc.ne.10) then
1097 C If input structure hasn't been supplied from the PDB file read or generate
1099 if (iranconf.eq.0 .and. .not. extconf) then
1100 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1101 & write (iout,'(a)') 'Initial geometry will be read in.'
1103 read(inp,'(8f10.5)',end=36,err=36)
1104 & ((c(l,k),l=1,3),k=1,nres),
1105 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1106 write (iout,*) "Exit READ_CART"
1107 write (iout,'(8f10.5)')
1108 & ((c(l,k),l=1,3),k=1,nres),
1109 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1110 call int_from_cart1(.true.)
1111 write (iout,*) "Finish INT_TO_CART"
1114 dc(j,i)=c(j,i+1)-c(j,i)
1115 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1119 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1121 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1122 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1128 write (iout,*) "Calling read_ang"
1129 call read_angles(inp,*36)
1130 write (iout,*) "Calling chainbuild"
1131 call chainbuild_extconf
1134 36 write (iout,'(a)') 'Error reading angle file.'
1136 call mpi_finalize( MPI_COMM_WORLD,IERR )
1138 stop 'Error reading angle file.'
1140 else if (extconf) then
1141 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1142 & write (iout,'(a)') 'Extended chain initial geometry.'
1144 theta(i)=90d0*deg2rad
1147 phi(i)=180d0*deg2rad
1150 alph(i)=110d0*deg2rad
1153 omeg(i)=-120d0*deg2rad
1154 if (itype(i).le.0) omeg(i)=-omeg(i)
1156 c from old chainbuild
1158 C Define the origin and orientation of the coordinate system and locate the
1159 C first three CA's and SC(2).
1163 * Build the alpha-carbon chain.
1166 call locate_next_res(i)
1169 C First and last SC must coincide with the corresponding CA.
1173 dc_norm(j,nres+1)=0.0D0
1174 dc(j,nres+nres)=0.0D0
1175 dc_norm(j,nres+nres)=0.0D0
1177 c(j,nres+nres)=c(j,nres)
1180 C Define the origin and orientation of the coordinate system and locate the
1181 C first three CA's and SC(2).
1185 * Build the alpha-carbon chain.
1188 call locate_next_res(i)
1191 C First and last SC must coincide with the corresponding CA.
1195 dc_norm(j,nres+1)=0.0D0
1196 dc(j,nres+nres)=0.0D0
1197 dc_norm(j,nres+nres)=0.0D0
1199 c(j,nres+nres)=c(j,nres)
1204 if(me.eq.king.or..not.out1file)
1205 & write (iout,'(a)') 'Random-generated initial geometry.'
1209 if (me.eq.king .or. fg_rank.eq.0 .and. (
1210 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1214 call gen_rand_conf(itmp,*30)
1216 30 write (iout,*) 'Failed to generate random conformation',
1217 & ', itrial=',itrial
1218 write (*,*) 'Processor:',me,
1219 & ' Failed to generate random conformation',
1229 write (iout,'(a,i3,a)') 'Processor:',me,
1230 & ' error in generating random conformation.'
1231 write (*,'(a,i3,a)') 'Processor:',me,
1232 & ' error in generating random conformation.'
1235 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1240 & ' error in generating random conformation.'
1245 elseif (modecalc.eq.4) then
1246 read (inp,'(a)') intinname
1247 open (intin,file=intinname,status='old',err=333)
1248 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1249 & write (iout,'(a)') 'intinname',intinname
1250 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1252 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1254 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1256 stop 'Error opening angle file.'
1260 C Generate distance constraints, if the PDB structure is to be regularized.
1261 if (nthread.gt.0) then
1262 call read_threadbase
1265 if (me.eq.king .or. .not. out1file)
1267 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1268 write (iout,'(/a,i3,a)')
1269 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1270 write (iout,'(20i4)') (iss(i),i=1,ns)
1272 write(iout,*)"Running with dynamic disulfide-bond formation"
1274 write (iout,'(/a/)') 'Pre-formed links are:'
1280 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1281 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1287 if (ns.gt.0.and.dyn_ss) then
1291 forcon(i-nss)=forcon(i)
1298 dyn_ss_mask(iss(i))=.true.
1301 if (i2ndstr.gt.0) call secstrp2dihc
1302 c call geom_to_var(nvar,x)
1303 c call etotal(energia(0))
1304 c call enerprint(energia(0))
1305 c call briefout(0,etot)
1307 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1308 cd write (iout,'(a)') 'Variable list:'
1309 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1311 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1312 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1313 & 'Processor',myrank,': end reading molecular data.'
1317 c--------------------------------------------------------------------------
1318 logical function seq_comp(itypea,itypeb,length)
1320 integer length,itypea(length),itypeb(length)
1323 if (itypea(i).ne.itypeb(i)) then
1331 c-----------------------------------------------------------------------------
1332 subroutine read_bridge
1333 C Read information about disulfide bridges.
1334 implicit real*8 (a-h,o-z)
1335 include 'DIMENSIONS'
1339 include 'COMMON.IOUNITS'
1340 include 'COMMON.GEO'
1341 include 'COMMON.VAR'
1342 include 'COMMON.INTERACT'
1343 include 'COMMON.LOCAL'
1344 include 'COMMON.NAMES'
1345 include 'COMMON.CHAIN'
1346 include 'COMMON.FFIELD'
1347 include 'COMMON.SBRIDGE'
1348 include 'COMMON.HEADER'
1349 include 'COMMON.CONTROL'
1350 include 'COMMON.DBASE'
1351 include 'COMMON.THREAD'
1352 include 'COMMON.TIME1'
1353 include 'COMMON.SETUP'
1354 C Read bridging residues.
1355 read (inp,*) ns,(iss(i),i=1,ns)
1357 if(me.eq.king.or..not.out1file)
1358 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1359 C Check whether the specified bridging residues are cystines.
1361 if (itype(iss(i)).ne.1) then
1362 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1363 & 'Do you REALLY think that the residue ',
1364 & restyp(itype(iss(i))),i,
1365 & ' can form a disulfide bridge?!!!'
1366 write (*,'(2a,i3,a)')
1367 & 'Do you REALLY think that the residue ',
1368 & restyp(itype(iss(i))),i,
1369 & ' can form a disulfide bridge?!!!'
1371 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1376 C Read preformed bridges.
1378 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1380 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1383 C Check if the residues involved in bridges are in the specified list of
1384 C bridging residues.
1387 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1388 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1389 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1390 & ' contains residues present in other pairs.'
1391 write (*,'(a,i3,a)') 'Disulfide pair',i,
1392 & ' contains residues present in other pairs.'
1394 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1400 if (ihpb(i).eq.iss(j)) goto 10
1402 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1405 if (jhpb(i).eq.iss(j)) goto 20
1407 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1413 ihpb(i)=ihpb(i)+nres
1414 jhpb(i)=jhpb(i)+nres
1420 c----------------------------------------------------------------------------
1421 subroutine read_x(kanal,*)
1422 implicit real*8 (a-h,o-z)
1423 include 'DIMENSIONS'
1424 include 'COMMON.GEO'
1425 include 'COMMON.VAR'
1426 include 'COMMON.CHAIN'
1427 include 'COMMON.IOUNITS'
1428 include 'COMMON.CONTROL'
1429 include 'COMMON.LOCAL'
1430 include 'COMMON.INTERACT'
1431 c Read coordinates from input
1433 read(kanal,'(8f10.5)',end=10,err=10)
1434 & ((c(l,k),l=1,3),k=1,nres),
1435 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1438 c(j,2*nres)=c(j,nres)
1440 call int_from_cart1(.false.)
1443 dc(j,i)=c(j,i+1)-c(j,i)
1444 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1448 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1450 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1451 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1459 c----------------------------------------------------------------------------
1460 subroutine read_threadbase
1461 implicit real*8 (a-h,o-z)
1462 include 'DIMENSIONS'
1463 include 'COMMON.IOUNITS'
1464 include 'COMMON.GEO'
1465 include 'COMMON.VAR'
1466 include 'COMMON.INTERACT'
1467 include 'COMMON.LOCAL'
1468 include 'COMMON.NAMES'
1469 include 'COMMON.CHAIN'
1470 include 'COMMON.FFIELD'
1471 include 'COMMON.SBRIDGE'
1472 include 'COMMON.HEADER'
1473 include 'COMMON.CONTROL'
1474 include 'COMMON.DBASE'
1475 include 'COMMON.THREAD'
1476 include 'COMMON.TIME1'
1477 C Read pattern database for threading.
1478 read (icbase,*) nseq
1480 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1481 & nres_base(2,i),nres_base(3,i)
1482 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1484 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1485 c & nres_base(2,i),nres_base(3,i)
1486 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1490 if (weidis.eq.0.0D0) weidis=0.1D0
1499 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1500 write (iout,'(a,i5)') 'nexcl: ',nexcl
1501 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1504 c------------------------------------------------------------------------------
1505 subroutine setup_var
1506 implicit real*8 (a-h,o-z)
1507 include 'DIMENSIONS'
1508 include 'COMMON.IOUNITS'
1509 include 'COMMON.GEO'
1510 include 'COMMON.VAR'
1511 include 'COMMON.INTERACT'
1512 include 'COMMON.LOCAL'
1513 include 'COMMON.NAMES'
1514 include 'COMMON.CHAIN'
1515 include 'COMMON.FFIELD'
1516 include 'COMMON.SBRIDGE'
1517 include 'COMMON.HEADER'
1518 include 'COMMON.CONTROL'
1519 include 'COMMON.DBASE'
1520 include 'COMMON.THREAD'
1521 include 'COMMON.TIME1'
1522 C Set up variable list.
1528 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1530 ialph(i,1)=nvar+nside
1534 if (indphi.gt.0) then
1536 else if (indback.gt.0) then
1541 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1544 c----------------------------------------------------------------------------
1545 subroutine gen_dist_constr
1546 C Generate CA distance constraints.
1547 implicit real*8 (a-h,o-z)
1548 include 'DIMENSIONS'
1549 include 'COMMON.IOUNITS'
1550 include 'COMMON.GEO'
1551 include 'COMMON.VAR'
1552 include 'COMMON.INTERACT'
1553 include 'COMMON.LOCAL'
1554 include 'COMMON.NAMES'
1555 include 'COMMON.CHAIN'
1556 include 'COMMON.FFIELD'
1557 include 'COMMON.SBRIDGE'
1558 include 'COMMON.HEADER'
1559 include 'COMMON.CONTROL'
1560 include 'COMMON.DBASE'
1561 include 'COMMON.THREAD'
1562 include 'COMMON.TIME1'
1563 dimension itype_pdb(maxres)
1564 common /pizda/ itype_pdb
1566 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1567 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1568 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1570 do i=nstart_sup,nstart_sup+nsup-1
1571 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1572 cd & ' seq_pdb', restyp(itype_pdb(i))
1573 do j=i+2,nstart_sup+nsup-1
1575 ihpb(nhpb)=i+nstart_seq-nstart_sup
1576 jhpb(nhpb)=j+nstart_seq-nstart_sup
1578 dhpb(nhpb)=dist(i,j)
1581 cd write (iout,'(a)') 'Distance constraints:'
1586 cd if (ii.gt.nres) then
1591 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1592 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1593 cd & dhpb(i),forcon(i)
1597 c----------------------------------------------------------------------------
1599 implicit real*8 (a-h,o-z)
1600 include 'DIMENSIONS'
1601 include 'COMMON.MAP'
1602 include 'COMMON.IOUNITS'
1603 character*3 angid(4) /'THE','PHI','ALP','OME'/
1604 character*80 mapcard,ucase
1606 read (inp,'(a)') mapcard
1607 mapcard=ucase(mapcard)
1608 if (index(mapcard,'PHI').gt.0) then
1610 else if (index(mapcard,'THE').gt.0) then
1612 else if (index(mapcard,'ALP').gt.0) then
1614 else if (index(mapcard,'OME').gt.0) then
1617 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1618 stop 'Error - illegal variable spec in MAP card.'
1620 call readi (mapcard,'RES1',res1(imap),0)
1621 call readi (mapcard,'RES2',res2(imap),0)
1622 if (res1(imap).eq.0) then
1623 res1(imap)=res2(imap)
1624 else if (res2(imap).eq.0) then
1625 res2(imap)=res1(imap)
1627 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1629 & 'Error - illegal definition of variable group in MAP.'
1630 stop 'Error - illegal definition of variable group in MAP.'
1632 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1633 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1634 call readi(mapcard,'NSTEP',nstep(imap),0)
1635 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1637 & 'Illegal boundary and/or step size specification in MAP.'
1638 stop 'Illegal boundary and/or step size specification in MAP.'
1643 c----------------------------------------------------------------------------
1645 implicit real*8 (a-h,o-z)
1646 include 'DIMENSIONS'
1647 include 'COMMON.IOUNITS'
1648 include 'COMMON.GEO'
1649 include 'COMMON.CSA'
1650 include 'COMMON.BANK'
1651 include 'COMMON.CONTROL'
1653 character*620 mcmcard
1654 call card_concat(mcmcard)
1656 call readi(mcmcard,'NCONF',nconf,50)
1657 call readi(mcmcard,'NADD',nadd,0)
1658 call readi(mcmcard,'JSTART',jstart,1)
1659 call readi(mcmcard,'JEND',jend,1)
1660 call readi(mcmcard,'NSTMAX',nstmax,500000)
1661 call readi(mcmcard,'N0',n0,1)
1662 call readi(mcmcard,'N1',n1,6)
1663 call readi(mcmcard,'N2',n2,4)
1664 call readi(mcmcard,'N3',n3,0)
1665 call readi(mcmcard,'N4',n4,0)
1666 call readi(mcmcard,'N5',n5,0)
1667 call readi(mcmcard,'N6',n6,10)
1668 call readi(mcmcard,'N7',n7,0)
1669 call readi(mcmcard,'N8',n8,0)
1670 call readi(mcmcard,'N9',n9,0)
1671 call readi(mcmcard,'N14',n14,0)
1672 call readi(mcmcard,'N15',n15,0)
1673 call readi(mcmcard,'N16',n16,0)
1674 call readi(mcmcard,'N17',n17,0)
1675 call readi(mcmcard,'N18',n18,0)
1677 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1679 call readi(mcmcard,'NDIFF',ndiff,2)
1680 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1681 call readi(mcmcard,'IS1',is1,1)
1682 call readi(mcmcard,'IS2',is2,8)
1683 call readi(mcmcard,'NRAN0',nran0,4)
1684 call readi(mcmcard,'NRAN1',nran1,2)
1685 call readi(mcmcard,'IRR',irr,1)
1686 call readi(mcmcard,'NSEED',nseed,20)
1687 call readi(mcmcard,'NTOTAL',ntotal,10000)
1688 call reada(mcmcard,'CUT1',cut1,2.0d0)
1689 call reada(mcmcard,'CUT2',cut2,5.0d0)
1690 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1691 call readi(mcmcard,'ICMAX',icmax,3)
1692 call readi(mcmcard,'IRESTART',irestart,0)
1693 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1696 call reada(mcmcard,'DELE',dele,20.0d0)
1697 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1698 call readi(mcmcard,'IREF',iref,0)
1699 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1700 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1701 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1702 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1703 write (iout,*) "NCONF_IN",nconf_in
1706 c----------------------------------------------------------------------------
1707 cfmc subroutine mcmfread
1708 cfmc implicit real*8 (a-h,o-z)
1709 cfmc include 'DIMENSIONS'
1710 cfmc include 'COMMON.MCMF'
1711 cfmc include 'COMMON.IOUNITS'
1712 cfmc include 'COMMON.GEO'
1713 cfmc character*80 ucase
1714 cfmc character*620 mcmcard
1715 cfmc call card_concat(mcmcard)
1717 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1718 cfmc write(iout,*)'MAXRANT=',maxrant
1719 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1720 cfmc write(iout,*)'MAXFAM=',maxfam
1721 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1722 cfmc write(iout,*)'NNET1=',nnet1
1723 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1724 cfmc write(iout,*)'NNET2=',nnet2
1725 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1726 cfmc write(iout,*)'NNET3=',nnet3
1727 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1728 cfmc write(iout,*)'ILASTT=',ilastt
1729 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1730 cfmc write(iout,*)'MAXSTR=',maxstr
1731 cfmc maxstr_f=maxstr/maxfam
1732 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1733 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1734 cfmc write(iout,*)'NMCMF=',nmcmf
1735 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1736 cfmc write(iout,*)'IFOCUS=',ifocus
1737 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1738 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1739 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1740 cfmc write(iout,*)'INTPRT=',intprt
1741 cfmc call readi(mcmcard,'IPRT',iprt,100)
1742 cfmc write(iout,*)'IPRT=',iprt
1743 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1744 cfmc write(iout,*)'IMAXTR=',imaxtr
1745 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1746 cfmc write(iout,*)'MAXEVEN=',maxeven
1747 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1748 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1749 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1750 cfmc write(iout,*)'INIMIN=',inimin
1751 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1752 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1753 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1754 cfmc write(iout,*)'NTHREAD=',nthread
1755 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1756 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1757 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1758 cfmc write(iout,*)'MAXPERT=',maxpert
1759 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1760 cfmc write(iout,*)'IRMSD=',irmsd
1761 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1762 cfmc write(iout,*)'DENEMIN=',denemin
1763 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1764 cfmc write(iout,*)'RCUT1S=',rcut1s
1765 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1766 cfmc write(iout,*)'RCUT1E=',rcut1e
1767 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1768 cfmc write(iout,*)'RCUT2S=',rcut2s
1769 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1770 cfmc write(iout,*)'RCUT2E=',rcut2e
1771 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1772 cfmc write(iout,*)'DPERT1=',d_pert1
1773 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1774 cfmc write(iout,*)'DPERT1A=',d_pert1a
1775 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1776 cfmc write(iout,*)'DPERT2=',d_pert2
1777 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1778 cfmc write(iout,*)'DPERT2A=',d_pert2a
1779 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1780 cfmc write(iout,*)'DPERT2B=',d_pert2b
1781 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1782 cfmc write(iout,*)'DPERT2C=',d_pert2c
1783 cfmc d_pert1=deg2rad*d_pert1
1784 cfmc d_pert1a=deg2rad*d_pert1a
1785 cfmc d_pert2=deg2rad*d_pert2
1786 cfmc d_pert2a=deg2rad*d_pert2a
1787 cfmc d_pert2b=deg2rad*d_pert2b
1788 cfmc d_pert2c=deg2rad*d_pert2c
1789 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1790 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1791 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1792 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1793 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1794 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1795 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1796 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1797 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1798 cfmc write(iout,*)'RCUTINI=',rcutini
1799 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1800 cfmc write(iout,*)'GRAT=',grat
1801 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1802 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1806 c----------------------------------------------------------------------------
1808 implicit real*8 (a-h,o-z)
1809 include 'DIMENSIONS'
1810 include 'COMMON.MCM'
1811 include 'COMMON.MCE'
1812 include 'COMMON.IOUNITS'
1814 character*320 mcmcard
1815 call card_concat(mcmcard)
1816 call readi(mcmcard,'MAXACC',maxacc,100)
1817 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1818 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1819 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1820 call readi(mcmcard,'MAXREPM',maxrepm,200)
1821 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1822 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1823 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1824 call reada(mcmcard,'E_UP',e_up,5.0D0)
1825 call reada(mcmcard,'DELTE',delte,0.1D0)
1826 call readi(mcmcard,'NSWEEP',nsweep,5)
1827 call readi(mcmcard,'NSTEPH',nsteph,0)
1828 call readi(mcmcard,'NSTEPC',nstepc,0)
1829 call reada(mcmcard,'TMIN',tmin,298.0D0)
1830 call reada(mcmcard,'TMAX',tmax,298.0D0)
1831 call readi(mcmcard,'NWINDOW',nwindow,0)
1832 call readi(mcmcard,'PRINT_MC',print_mc,0)
1833 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1834 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1835 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1836 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1837 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1838 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1839 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1840 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1841 if (nwindow.gt.0) then
1842 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1844 winlen(i)=winend(i)-winstart(i)+1
1847 if (tmax.lt.tmin) tmax=tmin
1848 if (tmax.eq.tmin) then
1852 if (nstepc.gt.0 .and. nsteph.gt.0) then
1853 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1854 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1856 C Probabilities of different move types
1857 sumpro_type(0)=0.0D0
1858 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1859 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1860 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1861 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1862 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1863 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1864 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1866 print *,'i',i,' sumprotype',sumpro_type(i)
1867 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1868 print *,'i',i,' sumprotype',sumpro_type(i)
1872 c----------------------------------------------------------------------------
1873 subroutine read_minim
1874 implicit real*8 (a-h,o-z)
1875 include 'DIMENSIONS'
1876 include 'COMMON.MINIM'
1877 include 'COMMON.IOUNITS'
1878 include 'COMMON.CONTROL'
1879 include 'COMMON.SETUP'
1881 character*320 minimcard
1882 call card_concat(minimcard)
1883 call readi(minimcard,'MAXMIN',maxmin,2000)
1884 call readi(minimcard,'MAXFUN',maxfun,5000)
1885 call readi(minimcard,'MINMIN',minmin,maxmin)
1886 call readi(minimcard,'MINFUN',minfun,maxmin)
1887 call reada(minimcard,'TOLF',tolf,1.0D-2)
1888 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1889 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1890 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1891 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1893 if (.not. out1file .or. me.eq.king) then
1895 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1896 & 'Options in energy minimization:'
1897 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1898 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1899 & 'MinMin:',MinMin,' MinFun:',MinFun,
1900 & ' TolF:',TolF,' RTolF:',RTolF
1906 c----------------------------------------------------------------------------
1907 subroutine read_angles(kanal,*)
1908 implicit real*8 (a-h,o-z)
1909 include 'DIMENSIONS'
1910 include 'COMMON.GEO'
1911 include 'COMMON.VAR'
1912 include 'COMMON.CHAIN'
1913 include 'COMMON.IOUNITS'
1914 include 'COMMON.CONTROL'
1915 c Read angles from input
1917 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1918 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1919 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1920 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1923 c 9/7/01 avoid 180 deg valence angle
1924 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1926 theta(i)=deg2rad*theta(i)
1927 phi(i)=deg2rad*phi(i)
1928 alph(i)=deg2rad*alph(i)
1929 omeg(i)=deg2rad*omeg(i)
1934 c----------------------------------------------------------------------------
1935 subroutine reada(rekord,lancuch,wartosc,default)
1937 character*(*) rekord,lancuch
1938 double precision wartosc,default
1941 iread=index(rekord,lancuch)
1942 if (iread.eq.0) then
1946 iread=iread+ilen(lancuch)+1
1947 read (rekord(iread:),*,err=10,end=10) wartosc
1952 c----------------------------------------------------------------------------
1953 subroutine readi(rekord,lancuch,wartosc,default)
1955 character*(*) rekord,lancuch
1956 integer wartosc,default
1959 iread=index(rekord,lancuch)
1960 if (iread.eq.0) then
1964 iread=iread+ilen(lancuch)+1
1965 read (rekord(iread:),*,err=10,end=10) wartosc
1970 c----------------------------------------------------------------------------
1971 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1974 integer tablica(dim),default
1975 character*(*) rekord,lancuch
1982 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1983 if (iread.eq.0) return
1984 iread=iread+ilen(lancuch)+1
1985 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1988 c----------------------------------------------------------------------------
1989 subroutine multreada(rekord,lancuch,tablica,dim,default)
1992 double precision tablica(dim),default
1993 character*(*) rekord,lancuch
2000 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2001 if (iread.eq.0) return
2002 iread=iread+ilen(lancuch)+1
2003 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2006 c----------------------------------------------------------------------------
2007 subroutine openunits
2008 implicit real*8 (a-h,o-z)
2009 include 'DIMENSIONS'
2012 character*16 form,nodename
2015 include 'COMMON.SETUP'
2016 include 'COMMON.IOUNITS'
2018 include 'COMMON.CONTROL'
2019 integer lenpre,lenpot,ilen,lentmp
2021 character*3 out1file_text,ucase
2024 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2025 call getenv_loc("PREFIX",prefix)
2027 call getenv_loc("POT",pot)
2028 call getenv_loc("DIRTMP",tmpdir)
2029 call getenv_loc("CURDIR",curdir)
2030 call getenv_loc("OUT1FILE",out1file_text)
2031 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2032 out1file_text=ucase(out1file_text)
2033 if (out1file_text(1:1).eq."Y") then
2036 out1file=fg_rank.gt.0
2041 if (lentmp.gt.0) then
2042 write (*,'(80(1h!))')
2043 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2044 write (*,'(80(1h!))')
2045 write (*,*)"All output files will be on node /tmp directory."
2047 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2048 if (me.eq.king) then
2049 write (*,*) "The master node is ",nodename
2050 else if (fg_rank.eq.0) then
2051 write (*,*) "I am the CG slave node ",nodename
2053 write (*,*) "I am the FG slave node ",nodename
2056 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2057 lenpre = lentmp+lenpre+1
2059 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2060 C Get the names and open the input files
2061 #if defined(WINIFL) || defined(WINPGI)
2062 open(1,file=pref_orig(:ilen(pref_orig))//
2063 & '.inp',status='old',readonly,shared)
2064 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2065 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2066 C Get parameter filenames and open the parameter files.
2067 call getenv_loc('BONDPAR',bondname)
2068 open (ibond,file=bondname,status='old',readonly,shared)
2069 call getenv_loc('THETPAR',thetname)
2070 open (ithep,file=thetname,status='old',readonly,shared)
2072 call getenv_loc('THETPARPDB',thetname_pdb)
2073 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2075 call getenv_loc('ROTPAR',rotname)
2076 open (irotam,file=rotname,status='old',readonly,shared)
2078 call getenv_loc('ROTPARPDB',rotname_pdb)
2079 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2081 call getenv_loc('TORPAR',torname)
2082 open (itorp,file=torname,status='old',readonly,shared)
2083 call getenv_loc('TORDPAR',tordname)
2084 open (itordp,file=tordname,status='old',readonly,shared)
2085 call getenv_loc('FOURIER',fouriername)
2086 open (ifourier,file=fouriername,status='old',readonly,shared)
2087 call getenv_loc('ELEPAR',elename)
2088 open (ielep,file=elename,status='old',readonly,shared)
2089 call getenv_loc('SIDEPAR',sidename)
2090 open (isidep,file=sidename,status='old',readonly,shared)
2091 #elif (defined CRAY) || (defined AIX)
2092 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2094 c print *,"Processor",myrank," opened file 1"
2095 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2096 c print *,"Processor",myrank," opened file 9"
2097 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2098 C Get parameter filenames and open the parameter files.
2099 call getenv_loc('BONDPAR',bondname)
2100 open (ibond,file=bondname,status='old',action='read')
2101 c print *,"Processor",myrank," opened file IBOND"
2102 call getenv_loc('THETPAR',thetname)
2103 open (ithep,file=thetname,status='old',action='read')
2104 c print *,"Processor",myrank," opened file ITHEP"
2106 call getenv_loc('THETPARPDB',thetname_pdb)
2107 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2109 call getenv_loc('ROTPAR',rotname)
2110 open (irotam,file=rotname,status='old',action='read')
2111 c print *,"Processor",myrank," opened file IROTAM"
2113 call getenv_loc('ROTPARPDB',rotname_pdb)
2114 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2116 call getenv_loc('TORPAR',torname)
2117 open (itorp,file=torname,status='old',action='read')
2118 c print *,"Processor",myrank," opened file ITORP"
2119 call getenv_loc('TORDPAR',tordname)
2120 open (itordp,file=tordname,status='old',action='read')
2121 c print *,"Processor",myrank," opened file ITORDP"
2122 call getenv_loc('SCCORPAR',sccorname)
2123 open (isccor,file=sccorname,status='old',action='read')
2124 c print *,"Processor",myrank," opened file ISCCOR"
2125 call getenv_loc('FOURIER',fouriername)
2126 open (ifourier,file=fouriername,status='old',action='read')
2127 c print *,"Processor",myrank," opened file IFOURIER"
2128 call getenv_loc('ELEPAR',elename)
2129 open (ielep,file=elename,status='old',action='read')
2130 c print *,"Processor",myrank," opened file IELEP"
2131 call getenv_loc('SIDEPAR',sidename)
2132 open (isidep,file=sidename,status='old',action='read')
2133 c print *,"Processor",myrank," opened file ISIDEP"
2134 c print *,"Processor",myrank," opened parameter files"
2136 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2137 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2138 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2139 C Get parameter filenames and open the parameter files.
2140 call getenv_loc('BONDPAR',bondname)
2141 open (ibond,file=bondname,status='old')
2142 call getenv_loc('THETPAR',thetname)
2143 open (ithep,file=thetname,status='old')
2145 call getenv_loc('THETPARPDB',thetname_pdb)
2146 open (ithep_pdb,file=thetname_pdb,status='old')
2148 call getenv_loc('ROTPAR',rotname)
2149 open (irotam,file=rotname,status='old')
2151 call getenv_loc('ROTPARPDB',rotname_pdb)
2152 open (irotam_pdb,file=rotname_pdb,status='old')
2154 call getenv_loc('TORPAR',torname)
2155 open (itorp,file=torname,status='old')
2156 call getenv_loc('TORDPAR',tordname)
2157 open (itordp,file=tordname,status='old')
2158 call getenv_loc('SCCORPAR',sccorname)
2159 open (isccor,file=sccorname,status='old')
2160 call getenv_loc('FOURIER',fouriername)
2161 open (ifourier,file=fouriername,status='old')
2162 call getenv_loc('ELEPAR',elename)
2163 open (ielep,file=elename,status='old')
2164 call getenv_loc('SIDEPAR',sidename)
2165 open (isidep,file=sidename,status='old')
2166 call getenv_loc('LIPTRANPAR',liptranname)
2167 open (iliptranpar,file=liptranname,status='old')
2169 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2171 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2172 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2173 C Get parameter filenames and open the parameter files.
2174 call getenv_loc('BONDPAR',bondname)
2175 open (ibond,file=bondname,status='old',readonly)
2176 call getenv_loc('THETPAR',thetname)
2177 open (ithep,file=thetname,status='old',readonly)
2178 call getenv_loc('ROTPAR',rotname)
2179 open (irotam,file=rotname,status='old',readonly)
2180 call getenv_loc('TORPAR',torname)
2181 open (itorp,file=torname,status='old',readonly)
2182 call getenv_loc('TORDPAR',tordname)
2183 open (itordp,file=tordname,status='old',readonly)
2184 call getenv_loc('SCCORPAR',sccorname)
2185 open (isccor,file=sccorname,status='old',readonly)
2187 call getenv_loc('THETPARPDB',thetname_pdb)
2188 c print *,"thetname_pdb ",thetname_pdb
2189 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2190 c print *,ithep_pdb," opened"
2192 call getenv_loc('FOURIER',fouriername)
2193 open (ifourier,file=fouriername,status='old',readonly)
2194 call getenv_loc('ELEPAR',elename)
2195 open (ielep,file=elename,status='old',readonly)
2196 call getenv_loc('SIDEPAR',sidename)
2197 open (isidep,file=sidename,status='old',readonly)
2198 call getenv_loc('LIPTRANPAR',liptranname)
2199 open (iliptranpar,file=liptranname,status='old',readonly)
2201 call getenv_loc('ROTPARPDB',rotname_pdb)
2202 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2207 C 8/9/01 In the newest version SCp interaction constants are read from a file
2208 C Use -DOLDSCP to use hard-coded constants instead.
2210 call getenv_loc('SCPPAR',scpname)
2211 #if defined(WINIFL) || defined(WINPGI)
2212 open (iscpp,file=scpname,status='old',readonly,shared)
2213 #elif (defined CRAY) || (defined AIX)
2214 open (iscpp,file=scpname,status='old',action='read')
2216 open (iscpp,file=scpname,status='old')
2218 open (iscpp,file=scpname,status='old',readonly)
2221 call getenv_loc('PATTERN',patname)
2222 #if defined(WINIFL) || defined(WINPGI)
2223 open (icbase,file=patname,status='old',readonly,shared)
2224 #elif (defined CRAY) || (defined AIX)
2225 open (icbase,file=patname,status='old',action='read')
2227 open (icbase,file=patname,status='old')
2229 open (icbase,file=patname,status='old',readonly)
2232 C Open output file only for CG processes
2233 c print *,"Processor",myrank," fg_rank",fg_rank
2234 if (fg_rank.eq.0) then
2236 if (nodes.eq.1) then
2239 npos = dlog10(dfloat(nodes-1))+1
2241 if (npos.lt.3) npos=3
2242 write (liczba,'(i1)') npos
2243 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2245 write (liczba,form) me
2246 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2247 & liczba(:ilen(liczba))
2248 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2250 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2252 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2253 & liczba(:ilen(liczba))//'.mol2'
2254 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2255 & liczba(:ilen(liczba))//'.stat'
2257 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2258 & //liczba(:ilen(liczba))//'.stat')
2259 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2262 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2263 & liczba(:ilen(liczba))//'.const'
2268 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2269 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2270 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2271 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2272 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2274 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2276 rest2name=prefix(:ilen(prefix))//'.rst'
2278 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2281 #if defined(AIX) || defined(PGI)
2282 if (me.eq.king .or. .not. out1file)
2283 & open(iout,file=outname,status='unknown')
2285 if (fg_rank.gt.0) then
2286 write (liczba,'(i3.3)') myrank/nfgtasks
2287 write (ll,'(bz,i3.3)') fg_rank
2288 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2293 open(igeom,file=intname,status='unknown',position='append')
2294 open(ipdb,file=pdbname,status='unknown')
2295 open(imol2,file=mol2name,status='unknown')
2296 open(istat,file=statname,status='unknown',position='append')
2298 c1out open(iout,file=outname,status='unknown')
2301 if (me.eq.king .or. .not.out1file)
2302 & open(iout,file=outname,status='unknown')
2304 if (fg_rank.gt.0) then
2305 write (liczba,'(i3.3)') myrank/nfgtasks
2306 write (ll,'(bz,i3.3)') fg_rank
2307 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2312 open(igeom,file=intname,status='unknown',access='append')
2313 open(ipdb,file=pdbname,status='unknown')
2314 open(imol2,file=mol2name,status='unknown')
2315 open(istat,file=statname,status='unknown',access='append')
2317 c1out open(iout,file=outname,status='unknown')
2320 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2321 csa_seed=prefix(:lenpre)//'.CSA.seed'
2322 csa_history=prefix(:lenpre)//'.CSA.history'
2323 csa_bank=prefix(:lenpre)//'.CSA.bank'
2324 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2325 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2326 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2327 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2328 csa_int=prefix(:lenpre)//'.int'
2329 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2330 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2331 csa_in=prefix(:lenpre)//'.CSA.in'
2332 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2335 write (iout,'(80(1h-))')
2336 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2337 write (iout,'(80(1h-))')
2338 write (iout,*) "Input file : ",
2339 & pref_orig(:ilen(pref_orig))//'.inp'
2340 write (iout,*) "Output file : ",
2341 & outname(:ilen(outname))
2343 write (iout,*) "Sidechain potential file : ",
2344 & sidename(:ilen(sidename))
2346 write (iout,*) "SCp potential file : ",
2347 & scpname(:ilen(scpname))
2349 write (iout,*) "Electrostatic potential file : ",
2350 & elename(:ilen(elename))
2351 write (iout,*) "Cumulant coefficient file : ",
2352 & fouriername(:ilen(fouriername))
2353 write (iout,*) "Torsional parameter file : ",
2354 & torname(:ilen(torname))
2355 write (iout,*) "Double torsional parameter file : ",
2356 & tordname(:ilen(tordname))
2357 write (iout,*) "SCCOR parameter file : ",
2358 & sccorname(:ilen(sccorname))
2359 write (iout,*) "Bond & inertia constant file : ",
2360 & bondname(:ilen(bondname))
2361 write (iout,*) "Bending parameter file : ",
2362 & thetname(:ilen(thetname))
2363 write (iout,*) "Rotamer parameter file : ",
2364 & rotname(:ilen(rotname))
2365 write (iout,*) "Thetpdb parameter file : ",
2366 & thetname_pdb(:ilen(thetname_pdb))
2367 write (iout,*) "Threading database : ",
2368 & patname(:ilen(patname))
2370 &write (iout,*)" DIRTMP : ",
2372 write (iout,'(80(1h-))')
2376 c----------------------------------------------------------------------------
2377 subroutine card_concat(card)
2378 implicit real*8 (a-h,o-z)
2379 include 'DIMENSIONS'
2380 include 'COMMON.IOUNITS'
2382 character*80 karta,ucase
2384 read (inp,'(a)') karta
2387 do while (karta(80:80).eq.'&')
2388 card=card(:ilen(card)+1)//karta(:79)
2389 read (inp,'(a)') karta
2392 card=card(:ilen(card)+1)//karta
2395 c----------------------------------------------------------------------------------
2397 implicit real*8 (a-h,o-z)
2398 include 'DIMENSIONS'
2399 include 'COMMON.CHAIN'
2400 include 'COMMON.IOUNITS'
2402 include 'COMMON.CONTROL'
2403 open(irest2,file=rest2name,status='unknown')
2404 read(irest2,*) totT,EK,potE,totE,t_bath
2407 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2410 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2412 if(usampl.or.homol_nset.gt.1) then
2413 read (irest2,*) iset
2418 c---------------------------------------------------------------------------------
2419 subroutine read_fragments
2420 implicit real*8 (a-h,o-z)
2421 include 'DIMENSIONS'
2425 include 'COMMON.SETUP'
2426 include 'COMMON.CHAIN'
2427 include 'COMMON.IOUNITS'
2429 include 'COMMON.CONTROL'
2431 read(inp,*) nset,nfrag,npair,nfrag_back
2432 if(me.eq.king.or..not.out1file)
2433 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2434 & " nfrag_back",nfrag_back
2436 read(inp,*) mset(iset1)
2438 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2440 if(me.eq.king.or..not.out1file)
2441 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2442 & ifrag(2,i,iset1), qinfrag(i,iset1)
2445 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2447 if(me.eq.king.or..not.out1file)
2448 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2449 & ipair(2,i,iset1), qinpair(i,iset1)
2452 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2453 & wfrag_back(3,i,iset1),
2454 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2455 if(me.eq.king.or..not.out1file)
2456 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2457 & wfrag_back(2,i,iset1),
2458 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2459 & ifrag_back(2,i,iset1)
2464 C---------------------------------------------------------------------------
2465 subroutine read_afminp
2466 implicit real*8 (a-h,o-z)
2467 include 'DIMENSIONS'
2471 include 'COMMON.SETUP'
2472 include 'COMMON.CONTROL'
2473 include 'COMMON.CHAIN'
2474 include 'COMMON.IOUNITS'
2475 include 'COMMON.SBRIDGE'
2476 character*320 afmcard
2477 c print *, "wchodze"
2478 call card_concat(afmcard)
2479 call readi(afmcard,"BEG",afmbeg,0)
2480 call readi(afmcard,"END",afmend,0)
2481 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2482 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2483 print *,'FORCE=' ,forceAFMconst
2484 CCCC NOW PROPERTIES FOR AFM
2487 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2489 distafminit=dsqrt(distafminit)
2490 print *,'initdist',distafminit
2494 c-------------------------------------------------------------------------------
2495 subroutine read_saxs_constr
2496 implicit real*8 (a-h,o-z)
2497 include 'DIMENSIONS'
2501 include 'COMMON.SETUP'
2502 include 'COMMON.CONTROL'
2503 include 'COMMON.CHAIN'
2504 include 'COMMON.IOUNITS'
2505 include 'COMMON.SBRIDGE'
2506 double precision cm(3)
2508 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2510 if (saxs_mode.eq.0) then
2511 c SAXS distance distribution
2513 read(inp,*) distsaxs(i),Psaxs(i)
2517 Cnorm = Cnorm + Psaxs(i)
2519 write (iout,*) "Cnorm",Cnorm
2521 Psaxs(i)=Psaxs(i)/Cnorm
2523 write (iout,*) "Normalized distance distribution from SAXS"
2525 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2529 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2531 write (iout,*) "Wsaxs0",Wsaxs0
2535 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2542 cm(j)=cm(j)+Csaxs(j,i)
2550 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2553 write (iout,*) "SAXS sphere coordinates"
2555 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2560 c-------------------------------------------------------------------------------
2561 subroutine read_dist_constr
2562 implicit real*8 (a-h,o-z)
2563 include 'DIMENSIONS'
2567 include 'COMMON.SETUP'
2568 include 'COMMON.CONTROL'
2569 include 'COMMON.CHAIN'
2570 include 'COMMON.IOUNITS'
2571 include 'COMMON.SBRIDGE'
2572 integer ifrag_(2,100),ipair_(2,1000)
2573 double precision wfrag_(100),wpair_(1000)
2574 character*5000 controlcard
2575 logical normalize,next
2577 double precision xlink(4,0:4) /
2579 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2580 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2581 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2582 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2583 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2584 c print *, "WCHODZE"
2585 write (iout,*) "Calling read_dist_constr"
2586 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2596 call card_concat(controlcard)
2597 next = index(controlcard,"NEXT").gt.0
2598 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2599 write (iout,*) "restr_type",restr_type
2600 call readi(controlcard,"NFRAG",nfrag_,0)
2601 call readi(controlcard,"NFRAG",nfrag_,0)
2602 call readi(controlcard,"NPAIR",npair_,0)
2603 call readi(controlcard,"NDIST",ndist_,0)
2604 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2605 if (restr_type.eq.10)
2606 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2607 if (restr_type.eq.12)
2608 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2609 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2610 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2611 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2612 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2613 normalize = index(controlcard,"NORMALIZE").gt.0
2614 write (iout,*) "WBOLTZD",wboltzd
2615 write (iout,*) "SCAL_PEAK",scal_peak
2616 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2617 write (iout,*) "IFRAG"
2619 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2621 write (iout,*) "IPAIR"
2623 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2625 if (nfrag_.gt.0) write (iout,*)
2626 & "Distance restraints as generated from reference structure"
2628 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2629 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2630 & ifrag_(2,i)=nstart_sup+nsup-1
2631 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2633 if (wfrag_(i).eq.0.0d0) cycle
2634 do j=ifrag_(1,i),ifrag_(2,i)-1
2635 do k=j+1,ifrag_(2,i)
2636 c write (iout,*) "j",j," k",k
2638 if (restr_type.eq.1) then
2644 forcon(nhpb)=wfrag_(i)
2645 else if (constr_dist.eq.2) then
2646 if (ddjk.le.dist_cut) then
2652 forcon(nhpb)=wfrag_(i)
2654 else if (restr_type.eq.3) then
2660 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2663 if (.not.out1file .or. me.eq.king)
2664 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2665 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2667 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2668 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2674 if (wpair_(i).eq.0.0d0) cycle
2682 do j=ifrag_(1,ii),ifrag_(2,ii)
2683 do k=ifrag_(1,jj),ifrag_(2,jj)
2685 if (restr_type.eq.1) then
2691 forcon(nhpb)=wpair_(i)
2692 else if (constr_dist.eq.2) then
2693 if (ddjk.le.dist_cut) then
2699 forcon(nhpb)=wpair_(i)
2701 else if (restr_type.eq.3) then
2707 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2710 if (.not.out1file .or. me.eq.king)
2711 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2712 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2714 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2715 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2722 write (iout,*) "Distance restraints as read from input"
2724 if (restr_type.eq.12) then
2725 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2726 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2727 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2728 & fordepth_peak(nhpb_peak+1),npeak
2729 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2730 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2731 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2732 c & fordepth_peak(nhpb_peak+1),npeak
2733 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2734 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2735 nhpb_peak=nhpb_peak+1
2736 irestr_type_peak(nhpb_peak)=12
2737 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2740 if (.not.out1file .or. me.eq.king)
2741 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2742 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2743 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2744 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2745 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2747 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2748 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2749 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2750 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2751 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2753 if (ibecarb_peak(nhpb_peak).gt.0) then
2754 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2755 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2757 else if (restr_type.eq.11) then
2758 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2759 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2760 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2761 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2763 irestr_type(nhpb)=11
2765 if (.not.out1file .or. me.eq.king)
2766 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2767 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2768 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2770 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2771 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2772 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2774 if (ibecarb(nhpb).gt.0) then
2775 ihpb(nhpb)=ihpb(nhpb)+nres
2776 jhpb(nhpb)=jhpb(nhpb)+nres
2778 else if (restr_type.eq.10) then
2779 c Cross-lonk Markov-like potential
2780 call card_concat(controlcard)
2781 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2782 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2784 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2785 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2786 if (index(controlcard,"ZL").gt.0) then
2788 else if (index(controlcard,"ADH").gt.0) then
2790 else if (index(controlcard,"PDH").gt.0) then
2792 else if (index(controlcard,"DSS").gt.0) then
2797 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2798 & xlink(1,link_type))
2799 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2800 & xlink(2,link_type))
2801 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2802 & xlink(3,link_type))
2803 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2804 & xlink(4,link_type))
2805 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2806 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2807 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2808 if (forcon(nhpb+1).le.0.0d0 .or.
2809 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2811 irestr_type(nhpb)=10
2812 if (ibecarb(nhpb).gt.0) then
2813 ihpb(nhpb)=ihpb(nhpb)+nres
2814 jhpb(nhpb)=jhpb(nhpb)+nres
2817 if (.not.out1file .or. me.eq.king)
2818 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2819 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2820 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2823 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2824 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2825 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2830 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2831 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2832 if (forcon(nhpb+1).gt.0.0d0) then
2834 if (dhpb1(nhpb).eq.0.0d0) then
2839 if (ibecarb(nhpb).gt.0) then
2840 ihpb(nhpb)=ihpb(nhpb)+nres
2841 jhpb(nhpb)=jhpb(nhpb)+nres
2843 if (dhpb(nhpb).eq.0.0d0)
2844 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2847 if (.not.out1file .or. me.eq.king)
2848 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2849 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2851 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2852 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2855 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2856 C if (forcon(nhpb+1).gt.0.0d0) then
2858 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2866 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2867 & fordepthmax=fordepth(i)
2870 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2875 c-------------------------------------------------------------------------------
2877 subroutine read_constr_homology
2879 include 'DIMENSIONS'
2883 include 'COMMON.SETUP'
2884 include 'COMMON.CONTROL'
2885 include 'COMMON.CHAIN'
2886 include 'COMMON.IOUNITS'
2888 include 'COMMON.GEO'
2889 include 'COMMON.INTERACT'
2890 include 'COMMON.NAMES'
2892 c For new homol impl
2894 include 'COMMON.VAR'
2897 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2899 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2900 c & sigma_odl_temp(maxres,maxres,max_template)
2902 character*24 model_ki_dist, model_ki_angle
2903 character*500 controlcard
2904 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2905 logical lprn /.true./
2910 c FP - Nov. 2014 Temporary specifications for new vars
2912 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2914 double precision, dimension (max_template,maxres) :: rescore
2915 double precision, dimension (max_template,maxres) :: rescore2
2916 double precision, dimension (max_template,maxres) :: rescore3
2917 character*24 pdbfile,tpl_k_rescore
2918 c -----------------------------------------------------------------
2919 c Reading multiple PDB ref structures and calculation of retraints
2920 c not using pre-computed ones stored in files model_ki_{dist,angle}
2922 c -----------------------------------------------------------------
2925 c Alternative: reading from input
2926 call card_concat(controlcard)
2927 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2928 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2929 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2930 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2931 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2932 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2933 dist1cut=(index(controlcard,'DIST1CUT').gt.0)
2934 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2935 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2936 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2937 if(.not.read2sigma.and.start_from_model) then
2938 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2939 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2940 start_from_model=.false.
2942 if(start_from_model .and. (me.eq.king .or. .not. out1file))
2943 & write(iout,*) 'START_FROM_MODELS is ON'
2944 if(start_from_model .and. rest) then
2945 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2946 write(iout,*) 'START_FROM_MODELS is OFF'
2947 write(iout,*) 'remove restart keyword from input'
2950 if (homol_nset.gt.1)then
2951 call card_concat(controlcard)
2952 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2953 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2954 write(iout,*) "iset homology_weight "
2956 write(iout,*) i,waga_homology(i)
2959 iset=mod(kolor,homol_nset)+1
2962 waga_homology(1)=1.0
2965 cd write (iout,*) "nnt",nnt," nct",nct
2972 c write(iout,*) 'nnt=',nnt,'nct=',nct
2975 do k=1,constr_homology
2988 if (read_homol_frag) then
2989 call read_klapaucjusz
2992 do k=1,constr_homology
2994 read(inp,'(a)') pdbfile
2995 if(me.eq.king .or. .not. out1file)
2996 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2997 & pdbfile(:ilen(pdbfile))
2998 open(ipdbin,file=pdbfile,status='old',err=33)
3000 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3001 & pdbfile(:ilen(pdbfile))
3004 c print *,'Begin reading pdb data'
3006 c Files containing res sim or local scores (former containing sigmas)
3009 write(kic2,'(bz,i2.2)') k
3011 tpl_k_rescore="template"//kic2//".sco"
3014 if (read2sigma) then
3015 call readpdb_template(k)
3020 c Distance restraints
3023 C Copy the coordinates from reference coordinates (?)
3027 c write (iout,*) "c(",j,i,") =",c(j,i)
3031 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3033 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3034 open (ientin,file=tpl_k_rescore,status='old')
3035 if (nnt.gt.1) rescore(k,1)=0.0d0
3036 do irec=nnt,nct ! loop for reading res sim
3037 if (read2sigma) then
3038 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3039 & rescore3_tmp,idomain_tmp
3041 idomain(k,i_tmp)=idomain_tmp
3042 rescore(k,i_tmp)=rescore_tmp
3043 rescore2(k,i_tmp)=rescore2_tmp
3044 rescore3(k,i_tmp)=rescore3_tmp
3045 if (.not. out1file .or. me.eq.king)
3046 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3047 & i_tmp,rescore2_tmp,rescore_tmp,
3048 & rescore3_tmp,idomain_tmp
3051 read (ientin,*,end=1401) rescore_tmp
3053 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3054 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3055 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3060 if (waga_dist.ne.0.0d0) then
3068 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3069 c write (iout,*) k,i,j,distal,dist2_cut
3071 if (dist1cut .and. k.gt.1) then
3073 if (l_homo(1,ii)) then
3079 sigma_odl(k,ii)=sigma_odl(1,ii)
3081 l_homo(k,ii)=.false.
3084 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3085 & .and. distal.le.dist2_cut ) then
3091 c write (iout,*) "k",k
3092 c write (iout,*) "i",i," j",j," constr_homology",
3097 if (read2sigma) then
3100 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3102 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3103 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3104 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3106 if (odl(k,ii).le.dist_cut) then
3107 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3110 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3111 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3113 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3114 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3118 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3121 l_homo(k,ii)=.false.
3129 c Theta, dihedral and SC retraints
3131 if (waga_angle.gt.0.0d0) then
3132 c open (ientin,file=tpl_k_sigma_dih,status='old')
3133 c do irec=1,maxres-3 ! loop for reading sigma_dih
3134 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3135 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3136 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3137 c & sigma_dih(k,i+nnt-1)
3142 if (idomain(k,i).eq.0) then
3146 dih(k,i)=phiref(i) ! right?
3147 c read (ientin,*) sigma_dih(k,i) ! original variant
3148 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3149 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3150 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3151 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3152 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3154 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3155 & rescore(k,i-2)+rescore(k,i-3))/4.0
3156 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3157 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3158 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3159 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3160 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3161 c Instead of res sim other local measure of b/b str reliability possible
3162 if (sigma_dih(k,i).ne.0)
3163 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3164 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3169 if (waga_theta.gt.0.0d0) then
3170 c open (ientin,file=tpl_k_sigma_theta,status='old')
3171 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3172 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3173 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3174 c & sigma_theta(k,i+nnt-1)
3179 do i = nnt+2,nct ! right? without parallel.
3180 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3181 c do i=ithet_start,ithet_end ! with FG parallel.
3182 if (idomain(k,i).eq.0) then
3183 sigma_theta(k,i)=0.0
3186 thetatpl(k,i)=thetaref(i)
3187 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3188 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3189 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3190 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3191 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3192 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3193 & rescore(k,i-2))/3.0
3194 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3195 if (sigma_theta(k,i).ne.0)
3196 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3198 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3199 c rescore(k,i-2) ! right expression ?
3200 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3204 if (waga_d.gt.0.0d0) then
3205 c open (ientin,file=tpl_k_sigma_d,status='old')
3206 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3207 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3208 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3209 c & sigma_d(k,i+nnt-1)
3213 do i = nnt,nct ! right? without parallel.
3214 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3215 c do i=loc_start,loc_end ! with FG parallel.
3216 if (itype(i).eq.10) cycle
3217 if (idomain(k,i).eq.0 ) then
3224 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3225 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3226 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3227 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3228 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3229 if (sigma_d(k,i).ne.0)
3230 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3232 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3233 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3234 c read (ientin,*) sigma_d(k,i) ! 1st variant
3239 c remove distance restraints not used in any model from the list
3240 c shift data in all arrays
3242 if (waga_dist.ne.0.0d0) then
3248 if (ii_in_use(ii).eq.0.and.liiflag) then
3252 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3253 & .not.liiflag.and.ii.eq.lim_odl) then
3254 if (ii.eq.lim_odl) then
3255 iishift=ii-iistart+1
3260 do ki=iistart,lim_odl-iishift
3261 ires_homo(ki)=ires_homo(ki+iishift)
3262 jres_homo(ki)=jres_homo(ki+iishift)
3263 ii_in_use(ki)=ii_in_use(ki+iishift)
3264 do k=1,constr_homology
3265 odl(k,ki)=odl(k,ki+iishift)
3266 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3267 l_homo(k,ki)=l_homo(k,ki+iishift)
3271 lim_odl=lim_odl-iishift
3277 endif ! .not. klapaucjusz
3279 if (constr_homology.gt.0) call homology_partition
3280 if (constr_homology.gt.0) call init_int_table
3281 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3282 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3286 if (.not.lprn) return
3287 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3288 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3289 write (iout,*) "Distance restraints from templates"
3291 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3292 & ii,ires_homo(ii),jres_homo(ii),
3293 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3294 & ki=1,constr_homology)
3296 write (iout,*) "Dihedral angle restraints from templates"
3298 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3299 & (rad2deg*dih(ki,i),
3300 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3302 write (iout,*) "Virtual-bond angle restraints from templates"
3304 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3305 & (rad2deg*thetatpl(ki,i),
3306 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3308 write (iout,*) "SC restraints from templates"
3310 write(iout,'(i5,100(4f8.2,4x))') i,
3311 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3312 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3315 c -----------------------------------------------------------------
3318 c----------------------------------------------------------------------
3321 subroutine flush(iu)
3326 subroutine flush(iu)
3331 c------------------------------------------------------------------------------
3332 subroutine copy_to_tmp(source)
3333 include "DIMENSIONS"
3334 include "COMMON.IOUNITS"
3335 character*(*) source
3336 character* 256 tmpfile
3340 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3341 inquire(file=tmpfile,exist=ex)
3343 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3344 & " to temporary directory..."
3345 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3346 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3350 c------------------------------------------------------------------------------
3351 subroutine move_from_tmp(source)
3352 include "DIMENSIONS"
3353 include "COMMON.IOUNITS"
3354 character*(*) source
3357 write (*,*) "Moving ",source(:ilen(source)),
3358 & " from temporary directory to working directory"
3359 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3360 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3363 c------------------------------------------------------------------------------
3364 subroutine random_init(seed)
3366 C Initialize random number generator
3368 implicit real*8 (a-h,o-z)
3369 include 'DIMENSIONS'
3375 logical OKRandom, prng_restart
3377 integer iseed_array(4)
3379 include 'COMMON.IOUNITS'
3380 include 'COMMON.TIME1'
3381 include 'COMMON.THREAD'
3382 include 'COMMON.SBRIDGE'
3383 include 'COMMON.CONTROL'
3384 include 'COMMON.MCM'
3385 include 'COMMON.MAP'
3386 include 'COMMON.HEADER'
3387 include 'COMMON.CSA'
3388 include 'COMMON.CHAIN'
3389 include 'COMMON.MUCA'
3391 include 'COMMON.FFIELD'
3392 include 'COMMON.SETUP'
3393 iseed=-dint(dabs(seed))
3394 if (iseed.eq.0) then
3395 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3396 & 'Random seed undefined. The program will stop.'
3397 write (*,'(/80(1h*)/20x,a/80(1h*))')
3398 & 'Random seed undefined. The program will stop.'
3400 call mpi_finalize(mpi_comm_world,ierr)
3402 stop 'Bad random seed.'
3405 if (fg_rank.eq.0) then
3409 if(me.eq.king .or. .not. out1file)
3410 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3411 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3412 OKRandom = prng_restart(me,iseedi8)
3415 tmp=65536.0d0**(4-i)
3416 iseed_array(i) = dint(seed/tmp)
3417 seed=seed-iseed_array(i)*tmp
3420 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3421 & (iseed_array(i),i=1,4)
3422 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3423 & (iseed_array(i),i=1,4)
3424 OKRandom = prng_restart(me,iseed_array)
3427 c r1 = prng_next(me)
3428 r1=ran_number(0.0D0,1.0D0)
3429 if(me.eq.king .or. .not. out1file)
3430 & write (iout,*) 'ran_num',r1
3431 if (r1.lt.0.0d0) OKRandom=.false.
3433 if (.not.OKRandom) then
3434 write (iout,*) 'PRNG IS NOT WORKING!!!'
3435 print *,'PRNG IS NOT WORKING!!!'
3438 call mpi_abort(mpi_comm_world,error_msg,ierr)
3441 write (iout,*) 'too many processors for parallel prng'
3442 write (*,*) 'too many processors for parallel prng'
3450 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3455 c -----------------------------------------------------------------
3456 subroutine read_klapaucjusz
3458 include 'DIMENSIONS'
3462 include 'COMMON.SETUP'
3463 include 'COMMON.CONTROL'
3464 include 'COMMON.CHAIN'
3465 include 'COMMON.IOUNITS'
3467 include 'COMMON.GEO'
3468 include 'COMMON.INTERACT'
3469 include 'COMMON.NAMES'
3470 character*256 fragfile
3471 integer ninclust(maxclust),inclust(max_template,maxclust),
3472 & nresclust(maxclust),iresclust(maxres,maxclust)
3475 character*24 model_ki_dist, model_ki_angle
3476 character*500 controlcard
3477 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3478 logical lprn /.true./
3484 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3485 double precision, dimension (max_template,maxres) :: rescore
3486 double precision, dimension (max_template,maxres) :: rescore2
3487 character*24 pdbfile,tpl_k_rescore
3490 c For new homol impl
3492 include 'COMMON.VAR'
3494 call getenv("FRAGFILE",fragfile)
3495 open(ientin,file=fragfile,status="old",err=10)
3496 read(ientin,*) constr_homology,nclust
3502 do k=1,constr_homology
3503 read(ientin,'(a)') pdbfile
3504 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3505 & pdbfile(:ilen(pdbfile))
3506 open(ipdbin,file=pdbfile,status='old',err=33)
3508 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3509 & pdbfile(:ilen(pdbfile))
3513 call readpdb_template(k)
3521 read(ientin,*) ninclust(i),nresclust(i)
3522 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3523 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3526 c Loop over clusters
3529 do ll = 1,ninclust(l)
3537 idomain(k,iresclust(i,l)+1) = 1
3539 idomain(k,iresclust(i,l)) = 1
3543 c Distance restraints
3546 C Copy the coordinates from reference coordinates (?)
3550 c write (iout,*) "c(",j,i,") =",c(j,i)
3553 call int_from_cart(.true.,.false.)
3554 call sc_loc_geom(.false.)
3556 thetaref(i)=theta(i)
3559 if (waga_dist.ne.0.0d0) then
3567 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3568 c write (iout,*) k,i,j,distal,dist2_cut
3570 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3571 & .and. distal.le.dist2_cut ) then
3577 c write (iout,*) "k",k
3578 c write (iout,*) "i",i," j",j," constr_homology",
3583 if (read2sigma) then
3586 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3588 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3589 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3590 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3592 if (odl(k,ii).le.dist_cut) then
3593 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3596 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3597 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3599 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3600 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3604 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3607 c l_homo(k,ii)=.false.
3614 c Theta, dihedral and SC retraints
3616 if (waga_angle.gt.0.0d0) then
3618 if (idomain(k,i).eq.0) then
3619 c sigma_dih(k,i)=0.0
3623 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3624 & rescore(k,i-2)+rescore(k,i-3))/4.0
3625 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3626 c & " sigma_dihed",sigma_dih(k,i)
3627 if (sigma_dih(k,i).ne.0)
3628 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3633 if (waga_theta.gt.0.0d0) then
3635 if (idomain(k,i).eq.0) then
3636 c sigma_theta(k,i)=0.0
3639 thetatpl(k,i)=thetaref(i)
3640 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3641 & rescore(k,i-2))/3.0
3642 if (sigma_theta(k,i).ne.0)
3643 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3647 if (waga_d.gt.0.0d0) then
3649 if (itype(i).eq.10) cycle
3650 if (idomain(k,i).eq.0 ) then
3657 sigma_d(k,i)=rescore(k,i)
3658 if (sigma_d(k,i).ne.0)
3659 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3660 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3666 c remove distance restraints not used in any model from the list
3667 c shift data in all arrays
3669 if (waga_dist.ne.0.0d0) then
3675 if (ii_in_use(ii).eq.0.and.liiflag) then
3679 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3680 & .not.liiflag.and.ii.eq.lim_odl) then
3681 if (ii.eq.lim_odl) then
3682 iishift=ii-iistart+1
3687 do ki=iistart,lim_odl-iishift
3688 ires_homo(ki)=ires_homo(ki+iishift)
3689 jres_homo(ki)=jres_homo(ki+iishift)
3690 ii_in_use(ki)=ii_in_use(ki+iishift)
3691 do k=1,constr_homology
3692 odl(k,ki)=odl(k,ki+iishift)
3693 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3694 l_homo(k,ki)=l_homo(k,ki+iishift)
3698 lim_odl=lim_odl-iishift
3705 10 stop "Error infragment file"
3707 c----------------------------------------------------------------------