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
116 read (INP,'(a)') titel
117 call card_concat(controlcard)
118 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
119 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
120 call reada(controlcard,'SEED',seed,0.0D0)
121 call random_init(seed)
122 C Set up the time limit (caution! The time must be input in minutes!)
123 read_cart=index(controlcard,'READ_CART').gt.0
124 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
125 write (iout,*) "constr_dist",constr_dist
126 call readi(controlcard,'NSAXS',nsaxs,0)
127 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
128 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
129 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
130 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
131 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
132 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
133 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
134 call readi(controlcard,'SYM',symetr,1)
135 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
136 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
137 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
138 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
139 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
140 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
141 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
142 call reada(controlcard,'DRMS',drms,0.1D0)
143 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
144 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
145 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
146 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
147 write (iout,'(a,f10.1)')'DRMS = ',drms
148 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
149 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
151 call readi(controlcard,'NZ_START',nz_start,0)
152 call readi(controlcard,'NZ_END',nz_end,0)
153 call readi(controlcard,'IZ_SC',iz_sc,0)
155 safety = 60.0d0*safety
158 call reada(controlcard,"T_BATH",t_bath,300.0d0)
159 minim=(index(controlcard,'MINIMIZE').gt.0)
160 dccart=(index(controlcard,'CART').gt.0)
161 overlapsc=(index(controlcard,'OVERLAP').gt.0)
162 overlapsc=.not.overlapsc
163 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
164 searchsc=.not.searchsc
165 sideadd=(index(controlcard,'SIDEADD').gt.0)
166 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
167 outpdb=(index(controlcard,'PDBOUT').gt.0)
168 outmol2=(index(controlcard,'MOL2OUT').gt.0)
169 pdbref=(index(controlcard,'PDBREF').gt.0)
170 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
171 indpdb=index(controlcard,'PDBSTART')
172 extconf=(index(controlcard,'EXTCONF').gt.0)
173 AFMlog=(index(controlcard,'AFM'))
174 selfguide=(index(controlcard,'SELFGUIDE'))
175 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
176 call readi(controlcard,'IPRINT',iprint,0)
177 call readi(controlcard,'MAXGEN',maxgen,10000)
178 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
179 call readi(controlcard,"KDIAG",kdiag,0)
180 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
181 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
182 & write (iout,*) "RESCALE_MODE",rescale_mode
183 split_ene=index(controlcard,'SPLIT_ENE').gt.0
184 if (index(controlcard,'REGULAR').gt.0.0D0) then
185 call reada(controlcard,'WEIDIS',weidis,0.1D0)
189 if (index(controlcard,'CHECKGRAD').gt.0) then
191 if (index(controlcard,'CART').gt.0) then
193 elseif (index(controlcard,'CARINT').gt.0) then
198 elseif (index(controlcard,'THREAD').gt.0) then
200 call readi(controlcard,'THREAD',nthread,0)
201 if (nthread.gt.0) then
202 call reada(controlcard,'WEIDIS',weidis,0.1D0)
205 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
206 stop 'Error termination in Read_Control.'
208 else if (index(controlcard,'MCMA').gt.0) then
210 else if (index(controlcard,'MCEE').gt.0) then
212 else if (index(controlcard,'MULTCONF').gt.0) then
214 else if (index(controlcard,'MAP').gt.0) then
216 call readi(controlcard,'MAP',nmap,0)
217 else if (index(controlcard,'CSA').gt.0) then
219 crc else if (index(controlcard,'ZSCORE').gt.0) then
221 crc ZSCORE is rm from UNRES, modecalc=9 is available
224 cfcm else if (index(controlcard,'MCMF').gt.0) then
226 else if (index(controlcard,'SOFTREG').gt.0) then
228 else if (index(controlcard,'CHECK_BOND').gt.0) then
230 else if (index(controlcard,'TEST').gt.0) then
232 else if (index(controlcard,'MD').gt.0) then
234 else if (index(controlcard,'RE ').gt.0) then
238 lmuca=index(controlcard,'MUCA').gt.0
239 call readi(controlcard,'MUCADYN',mucadyn,0)
240 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
241 if (lmuca .and. (me.eq.king .or. .not.out1file ))
243 write (iout,*) 'MUCADYN=',mucadyn
244 write (iout,*) 'MUCASMOOTH=',muca_smooth
247 iscode=index(controlcard,'ONE_LETTER')
248 indphi=index(controlcard,'PHI')
249 indback=index(controlcard,'BACK')
250 iranconf=index(controlcard,'RAND_CONF')
251 i2ndstr=index(controlcard,'USE_SEC_PRED')
252 gradout=index(controlcard,'GRADOUT').gt.0
253 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
254 C DISTCHAINMAX become obsolete for periodic boundry condition
255 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
256 C Reading the dimensions of box in x,y,z coordinates
257 call reada(controlcard,'BOXX',boxxsize,100.0d0)
258 call reada(controlcard,'BOXY',boxysize,100.0d0)
259 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
260 c Cutoff range for interactions
261 call reada(controlcard,"R_CUT",r_cut,15.0d0)
262 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
263 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
264 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
265 if (lipthick.gt.0.0d0) then
266 bordliptop=(boxzsize+lipthick)/2.0
267 bordlipbot=bordliptop-lipthick
269 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
270 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
271 buflipbot=bordlipbot+lipbufthick
272 bufliptop=bordliptop-lipbufthick
273 if ((lipbufthick*2.0d0).gt.lipthick)
274 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
276 c write(iout,*) "bordliptop=",bordliptop
277 c write(iout,*) "bordlipbot=",bordlipbot
278 c write(iout,*) "bufliptop=",bufliptop
279 c write(iout,*) "buflipbot=",buflipbot
282 if (me.eq.king .or. .not.out1file )
283 & write (iout,*) "DISTCHAINMAX",distchainmax
285 if(me.eq.king.or..not.out1file)
286 & write (iout,'(2a)') diagmeth(kdiag),
287 & ' routine used to diagonalize matrices.'
290 c--------------------------------------------------------------------------
291 subroutine read_REMDpar
295 implicit real*8 (a-h,o-z)
297 include 'COMMON.IOUNITS'
298 include 'COMMON.TIME1'
301 include 'COMMON.LANGEVIN'
303 include 'COMMON.LANGEVIN.lang0'
305 include 'COMMON.INTERACT'
306 include 'COMMON.NAMES'
308 include 'COMMON.REMD'
309 include 'COMMON.CONTROL'
310 include 'COMMON.SETUP'
312 character*320 controlcard
313 character*3200 controlcard1
314 integer iremd_m_total
316 if(me.eq.king.or..not.out1file)
317 & write (iout,*) "REMD setup"
319 call card_concat(controlcard)
320 call readi(controlcard,"NREP",nrep,3)
321 call readi(controlcard,"NSTEX",nstex,1000)
322 call reada(controlcard,"RETMIN",retmin,10.0d0)
323 call reada(controlcard,"RETMAX",retmax,1000.0d0)
324 mremdsync=(index(controlcard,'SYNC').gt.0)
325 call readi(controlcard,"NSYN",i_sync_step,100)
326 restart1file=(index(controlcard,'REST1FILE').gt.0)
327 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
328 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
329 if(max_cache_traj_use.gt.max_cache_traj)
330 & max_cache_traj_use=max_cache_traj
331 if(me.eq.king.or..not.out1file) then
332 cd if (traj1file) then
333 crc caching is in testing - NTWX is not ignored
334 cd write (iout,*) "NTWX value is ignored"
335 cd write (iout,*) " trajectory is stored to one file by master"
336 cd write (iout,*) " before exchange at NSTEX intervals"
338 write (iout,*) "NREP= ",nrep
339 write (iout,*) "NSTEX= ",nstex
340 write (iout,*) "SYNC= ",mremdsync
341 write (iout,*) "NSYN= ",i_sync_step
342 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
345 if (index(controlcard,'TLIST').gt.0) then
347 call card_concat(controlcard1)
348 read(controlcard1,*) (remd_t(i),i=1,nrep)
349 if(me.eq.king.or..not.out1file)
350 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
353 if (index(controlcard,'MLIST').gt.0) then
355 call card_concat(controlcard1)
356 read(controlcard1,*) (remd_m(i),i=1,nrep)
357 if(me.eq.king.or..not.out1file) then
358 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
361 iremd_m_total=iremd_m_total+remd_m(i)
363 write (iout,*) 'Total number of replicas ',iremd_m_total
366 if(me.eq.king.or..not.out1file)
367 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
370 c--------------------------------------------------------------------------
371 subroutine read_MDpar
375 implicit real*8 (a-h,o-z)
377 include 'COMMON.IOUNITS'
378 include 'COMMON.TIME1'
381 include 'COMMON.LANGEVIN'
383 include 'COMMON.LANGEVIN.lang0'
385 include 'COMMON.INTERACT'
386 include 'COMMON.NAMES'
388 include 'COMMON.SETUP'
389 include 'COMMON.CONTROL'
390 include 'COMMON.SPLITELE'
392 character*320 controlcard
394 call card_concat(controlcard)
395 call readi(controlcard,"NSTEP",n_timestep,1000000)
396 call readi(controlcard,"NTWE",ntwe,100)
397 call readi(controlcard,"NTWX",ntwx,1000)
398 call reada(controlcard,"DT",d_time,1.0d-1)
399 call reada(controlcard,"DVMAX",dvmax,2.0d1)
400 call reada(controlcard,"DAMAX",damax,1.0d1)
401 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
402 call readi(controlcard,"LANG",lang,0)
403 RESPA = index(controlcard,"RESPA") .gt. 0
404 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
405 ntime_split0=ntime_split
406 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
407 ntime_split0=ntime_split
408 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
409 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
410 rest = index(controlcard,"REST").gt.0
411 tbf = index(controlcard,"TBF").gt.0
412 usampl = index(controlcard,"USAMPL").gt.0
414 mdpdb = index(controlcard,"MDPDB").gt.0
415 call reada(controlcard,"T_BATH",t_bath,300.0d0)
416 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
417 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
418 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
419 if (count_reset_moment.eq.0) count_reset_moment=1000000000
420 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
421 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
422 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
423 if (count_reset_vel.eq.0) count_reset_vel=1000000000
424 large = index(controlcard,"LARGE").gt.0
425 print_compon = index(controlcard,"PRINT_COMPON").gt.0
426 rattle = index(controlcard,"RATTLE").gt.0
427 preminim = index(controlcard,"PREMINIM").gt.0
429 dccart=(index(controlcard,'CART').gt.0)
432 c if performing umbrella sampling, fragments constrained are read from the fragment file
438 if(me.eq.king.or..not.out1file) then
440 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
442 write (iout,'(a)') "The units are:"
443 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
444 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
445 & " acceleration: angstrom/(48.9 fs)**2"
446 write (iout,'(a)') "energy: kcal/mol, temperature: K"
448 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
449 write (iout,'(a60,f10.5,a)')
450 & "Initial time step of numerical integration:",d_time,
452 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
454 write (iout,'(2a,i4,a)')
455 & "A-MTS algorithm used; initial time step for fast-varying",
456 & " short-range forces split into",ntime_split," steps."
457 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
458 & r_cut," lambda",rlamb
460 write (iout,'(2a,f10.5)')
461 & "Maximum acceleration threshold to reduce the time step",
462 & "/increase split number:",damax
463 write (iout,'(2a,f10.5)')
464 & "Maximum predicted energy drift to reduce the timestep",
465 & "/increase split number:",edriftmax
466 write (iout,'(a60,f10.5)')
467 & "Maximum velocity threshold to reduce velocities:",dvmax
468 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
469 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
470 if (rattle) write (iout,'(a60)')
471 & "Rattle algorithm used to constrain the virtual bonds"
472 if (preminim .or. iranconf.gt.0) then
474 & "Initial structure will be energy-minimized"
479 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
480 call reada(controlcard,"RWAT",rwat,1.4d0)
481 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
482 surfarea=index(controlcard,"SURFAREA").gt.0
483 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
484 if(me.eq.king.or..not.out1file)then
485 write (iout,'(/a,$)') "Langevin dynamics calculation"
488 & " with direct integration of Langevin equations"
489 else if (lang.eq.2) then
490 write (iout,'(a/)') " with TINKER stochasic MD integrator"
491 else if (lang.eq.3) then
492 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
493 else if (lang.eq.4) then
494 write (iout,'(a/)') " in overdamped mode"
496 write (iout,'(//a,i5)')
497 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
500 write (iout,'(a60,f10.5)') "Temperature:",t_bath
501 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
502 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
503 write (iout,'(a60,f10.5)')
504 & "Scaling factor of the friction forces:",scal_fric
505 if (surfarea) write (iout,'(2a,i10,a)')
506 & "Friction coefficients will be scaled by solvent-accessible",
507 & " surface area every",reset_fricmat," steps."
509 c Calculate friction coefficients and bounds of stochastic forces
510 eta=6*pi*cPoise*etawat
511 if(me.eq.king.or..not.out1file)
512 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
514 gamp=scal_fric*(pstok+rwat)*eta
515 stdfp=dsqrt(2*Rb*t_bath/d_time)
517 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
518 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
520 if(me.eq.king.or..not.out1file)then
521 write (iout,'(/2a/)')
522 & "Radii of site types and friction coefficients and std's of",
523 & " stochastic forces of fully exposed sites"
524 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
526 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
527 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
531 if(me.eq.king.or..not.out1file)then
532 write (iout,'(a)') "Berendsen bath calculation"
533 write (iout,'(a60,f10.5)') "Temperature:",t_bath
534 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
536 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
537 & count_reset_moment," steps"
539 & write (iout,'(a,i10,a)')
540 & "Velocities will be reset at random every",count_reset_vel,
544 if(me.eq.king.or..not.out1file)
545 & write (iout,'(a31)') "Microcanonical mode calculation"
547 if(me.eq.king.or..not.out1file)then
548 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
550 write(iout,*) "MD running with constraints."
551 write(iout,*) "Equilibration time ", eq_time, " mtus."
552 write(iout,*) "Constraining ", nfrag," fragments."
553 write(iout,*) "Length of each fragment, weight and q0:"
555 write (iout,*) "Set of restraints #",iset
557 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
558 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
560 write(iout,*) "constraints between ", npair, "fragments."
561 write(iout,*) "constraint pairs, weights and q0:"
563 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
564 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
566 write(iout,*) "angle constraints within ", nfrag_back,
567 & "backbone fragments."
568 write(iout,*) "fragment, weights:"
570 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
571 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
572 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
575 iset=mod(kolor,nset)+1
578 if(me.eq.king.or..not.out1file)
579 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
582 c------------------------------------------------------------------------------
585 C Read molecular data.
587 implicit real*8 (a-h,o-z)
593 include 'COMMON.IOUNITS'
596 include 'COMMON.INTERACT'
597 include 'COMMON.LOCAL'
598 include 'COMMON.NAMES'
599 include 'COMMON.CHAIN'
600 include 'COMMON.FFIELD'
601 include 'COMMON.SBRIDGE'
602 include 'COMMON.HEADER'
603 include 'COMMON.CONTROL'
604 include 'COMMON.DBASE'
605 include 'COMMON.THREAD'
606 include 'COMMON.CONTACTS'
607 include 'COMMON.TORCNSTR'
608 include 'COMMON.TIME1'
609 include 'COMMON.BOUNDS'
611 include 'COMMON.SETUP'
612 character*4 sequence(maxres)
614 double precision x(maxvar)
615 character*256 pdbfile
616 character*320 weightcard
617 character*80 weightcard_t,ucase
618 dimension itype_pdb(maxres)
619 common /pizda/ itype_pdb
620 logical seq_comp,fail
621 double precision energia(0:n_ene)
627 C Read weights of the subsequent energy terms.
628 call card_concat(weightcard)
629 call reada(weightcard,'WLONG',wlong,1.0D0)
630 call reada(weightcard,'WSC',wsc,wlong)
631 call reada(weightcard,'WSCP',wscp,wlong)
632 call reada(weightcard,'WELEC',welec,1.0D0)
633 call reada(weightcard,'WVDWPP',wvdwpp,welec)
634 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
635 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
636 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
637 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
638 call reada(weightcard,'WTURN3',wturn3,1.0D0)
639 call reada(weightcard,'WTURN4',wturn4,1.0D0)
640 call reada(weightcard,'WTURN6',wturn6,1.0D0)
641 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
642 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
643 call reada(weightcard,'WBOND',wbond,1.0D0)
644 call reada(weightcard,'WTOR',wtor,1.0D0)
645 call reada(weightcard,'WTORD',wtor_d,1.0D0)
646 call reada(weightcard,'WANG',wang,1.0D0)
647 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
648 call reada(weightcard,'SCAL14',scal14,0.4D0)
649 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
650 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
651 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
652 call reada(weightcard,'TEMP0',temp0,300.0d0)
653 call reada(weightcard,'WLT',wliptran,0.0D0)
654 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
655 if (index(weightcard,'SOFT').gt.0) ipot=6
656 C 12/1/95 Added weight for the multi-body term WCORR
657 call reada(weightcard,'WCORRH',wcorr,1.0D0)
658 if (wcorr4.gt.0.0d0) wcorr=wcorr4
679 if(me.eq.king.or..not.out1file)
680 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
681 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
683 10 format (/'Energy-term weights (unscaled):'//
684 & 'WSCC= ',f10.6,' (SC-SC)'/
685 & 'WSCP= ',f10.6,' (SC-p)'/
686 & 'WELEC= ',f10.6,' (p-p electr)'/
687 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
688 & 'WBOND= ',f10.6,' (stretching)'/
689 & 'WANG= ',f10.6,' (bending)'/
690 & 'WSCLOC= ',f10.6,' (SC local)'/
691 & 'WTOR= ',f10.6,' (torsional)'/
692 & 'WTORD= ',f10.6,' (double torsional)'/
693 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
694 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
695 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
696 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
697 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
698 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
699 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
700 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
701 & 'WTURN6= ',f10.6,' (turns, 6th order)')
702 if(me.eq.king.or..not.out1file)then
703 if (wcorr4.gt.0.0d0) then
704 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
705 & 'between contact pairs of peptide groups'
706 write (iout,'(2(a,f5.3/))')
707 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
708 & 'Range of quenching the correlation terms:',2*delt_corr
709 else if (wcorr.gt.0.0d0) then
710 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
711 & 'between contact pairs of peptide groups'
713 write (iout,'(a,f8.3)')
714 & 'Scaling factor of 1,4 SC-p interactions:',scal14
715 write (iout,'(a,f8.3)')
716 & 'General scaling factor of SC-p interactions:',scalscp
718 r0_corr=cutoff_corr-delt_corr
720 aad(i,1)=scalscp*aad(i,1)
721 aad(i,2)=scalscp*aad(i,2)
722 bad(i,1)=scalscp*bad(i,1)
723 bad(i,2)=scalscp*bad(i,2)
725 call rescale_weights(t_bath)
726 if(me.eq.king.or..not.out1file)
727 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
728 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
730 22 format (/'Energy-term weights (scaled):'//
731 & 'WSCC= ',f10.6,' (SC-SC)'/
732 & 'WSCP= ',f10.6,' (SC-p)'/
733 & 'WELEC= ',f10.6,' (p-p electr)'/
734 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
735 & 'WBOND= ',f10.6,' (stretching)'/
736 & 'WANG= ',f10.6,' (bending)'/
737 & 'WSCLOC= ',f10.6,' (SC local)'/
738 & 'WTOR= ',f10.6,' (torsional)'/
739 & 'WTORD= ',f10.6,' (double torsional)'/
740 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
741 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
742 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
743 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
744 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
745 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
746 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
747 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
748 & 'WTURN6= ',f10.6,' (turns, 6th order)')
749 if(me.eq.king.or..not.out1file)
750 & write (iout,*) "Reference temperature for weights calculation:",
752 call reada(weightcard,"D0CM",d0cm,3.78d0)
753 call reada(weightcard,"AKCM",akcm,15.1d0)
754 call reada(weightcard,"AKTH",akth,11.0d0)
755 call reada(weightcard,"AKCT",akct,12.0d0)
756 call reada(weightcard,"V1SS",v1ss,-1.08d0)
757 call reada(weightcard,"V2SS",v2ss,7.61d0)
758 call reada(weightcard,"V3SS",v3ss,13.7d0)
759 call reada(weightcard,"EBR",ebr,-5.50D0)
760 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
762 dyn_ss_mask(i)=.false.
766 dyn_ssbond_ij(i,j)=1.0d300
769 call reada(weightcard,"HT",Ht,0.0D0)
771 ss_depth=ebr/wsc-0.25*eps(1,1)
772 Ht=Ht/wsc-0.25*eps(1,1)
773 akcm=akcm*wstrain/wsc
774 akth=akth*wstrain/wsc
775 akct=akct*wstrain/wsc
776 v1ss=v1ss*wstrain/wsc
777 v2ss=v2ss*wstrain/wsc
778 v3ss=v3ss*wstrain/wsc
780 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
783 if(me.eq.king.or..not.out1file) then
784 write (iout,*) "Parameters of the SS-bond potential:"
785 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
787 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
788 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
789 write (iout,*)" HT",Ht
790 print *,'indpdb=',indpdb,' pdbref=',pdbref
792 if (indpdb.gt.0 .or. pdbref) then
793 read(inp,'(a)') pdbfile
794 if(me.eq.king.or..not.out1file)
795 & write (iout,'(2a)') 'PDB data will be read from file ',
796 & pdbfile(:ilen(pdbfile))
797 open(ipdbin,file=pdbfile,status='old',err=33)
799 33 write (iout,'(a)') 'Error opening PDB file.'
802 c print *,'Begin reading pdb data'
811 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
812 & (crefjlee(j,i+nres),j=1,3)
815 c print *,'Finished reading pdb data'
816 if(me.eq.king.or..not.out1file)
817 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
818 & ' nstart_sup=',nstart_sup
820 itype_pdb(i)=itype(i)
824 nct=nstart_sup+nsup-1
825 call contact(.false.,ncont_ref,icont_ref,co)
828 if(me.eq.king.or..not.out1file)
829 & write(iout,*)'Adding sidechains'
833 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
836 do while (fail.and.nsi.le.maxsi)
837 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
841 c Calculalte only the coordinates of the current sidechain; no need to rebuild
843 call locate_side_chain(i)
844 if(fail) write(iout,*)'Adding sidechain failed for res ',
845 & i,' after ',nsi,' trials'
848 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
852 if (indpdb.eq.0) then
853 C Read sequence if not taken from the pdb file.
855 c print *,'nres=',nres
856 if (iscode.gt.0) then
857 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
859 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
861 C Convert sequence to numeric code
863 itype(i)=rescode(i,sequence(i),iscode)
865 C Assign initial virtual bond lengths
871 vbld(i+nres)=dsc(iabs(itype(i)))
872 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
873 c write (iout,*) "i",i," itype",itype(i),
874 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
878 c print '(20i4)',(itype(i),i=1,nres)
881 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
883 if (itype(i).eq.ntyp1) then
887 else if (iabs(itype(i+1)).ne.20) then
889 else if (iabs(itype(i)).ne.20) then
896 if(me.eq.king.or..not.out1file)then
897 write (iout,*) "ITEL"
899 write (iout,*) i,itype(i),itel(i)
901 print *,'Call Read_Bridge.'
904 C 8/13/98 Set limits to generating the dihedral angles
909 read (inp,*) ndih_constr
910 write (iout,*) "ndish_constr",ndih_constr
911 if (ndih_constr.gt.0) then
913 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
914 if(me.eq.king.or..not.out1file)then
916 & 'There are',ndih_constr,' constraints on phi angles.'
918 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
922 phi0(i)=deg2rad*phi0(i)
923 drange(i)=deg2rad*drange(i)
925 if(me.eq.king.or..not.out1file)
926 & write (iout,*) 'FTORS',ftors
929 phibound(1,ii) = phi0(i)-drange(i)
930 phibound(2,ii) = phi0(i)+drange(i)
937 write (iout,'(a)') 'Boundaries in phi angle sampling:'
939 write (iout,'(a3,i5,2f10.1)')
940 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
946 cd print *,'NNT=',NNT,' NCT=',NCT
947 if (itype(1).eq.ntyp1) nnt=2
948 if (itype(nres).eq.ntyp1) nct=nct-1
950 if(me.eq.king.or..not.out1file)
951 & write (iout,'(a,i3)') 'nsup=',nsup
953 if (nsup.le.(nct-nnt+1)) then
954 do i=0,nct-nnt+1-nsup
955 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
961 & 'Error - sequences to be superposed do not match.'
964 do i=0,nsup-(nct-nnt+1)
965 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
967 nstart_sup=nstart_sup+i
973 & 'Error - sequences to be superposed do not match.'
976 if (nsup.eq.0) nsup=nct-nnt
977 if (nstart_sup.eq.0) nstart_sup=nnt
978 if (nstart_seq.eq.0) nstart_seq=nnt
979 if(me.eq.king.or..not.out1file)
980 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
981 & ' nstart_seq=',nstart_seq
983 c--- Zscore rms -------
984 if (nz_start.eq.0) nz_start=nnt
985 if (nz_end.eq.0 .and. nsup.gt.0) then
987 else if (nz_end.eq.0) then
990 if(me.eq.king.or..not.out1file)then
991 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
992 write (iout,*) 'IZ_SC=',iz_sc
994 c----------------------
997 if (.not.pdbref) then
998 call read_angles(inp,*38)
1000 38 write (iout,'(a)') 'Error reading reference structure.'
1002 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1003 stop 'Error reading reference structure'
1005 39 call chainbuild_extconf
1007 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1014 cref(j,i,kkk)=c(j,i)
1017 call contact(.true.,ncont_ref,icont_ref,co)
1021 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1023 if (constr_dist.gt.0) call read_dist_constr
1024 write (iout,*) "After read_dist_constr nhpb",nhpb
1025 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1026 if(me.eq.king.or..not.out1file)
1027 & write (iout,*) 'Contact order:',co
1029 if(me.eq.king.or..not.out1file)
1030 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1033 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1035 if(me.eq.king.or..not.out1file)
1036 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1037 & icont_ref(1,i),' ',
1038 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1041 write (iout,*) "calling read_saxs_consrtr",nsaxs
1042 if (nsaxs.gt.0) call read_saxs_constr
1044 if (constr_homology.gt.0) then
1045 call read_constr_homology
1046 if (indpdb.gt.0 .or. pdbref) then
1049 c(j,i)=crefjlee(j,i)
1050 cref(j,i,1)=crefjlee(j,i)
1055 write (iout,*) "Array C"
1057 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1058 & (c(j,i+nres),j=1,3)
1060 write (iout,*) "Array Cref"
1062 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
1063 & (cref(j,i+nres,1),j=1,3)
1066 call int_from_cart1(.false.)
1067 call sc_loc_geom(.false.)
1069 thetaref(i)=theta(i)
1074 dc(j,i)=c(j,i+1)-c(j,i)
1075 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1080 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1081 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1089 if (nhpb.gt.0) call hpb_partition
1090 if (peak.gt.0) call NMRpeak_partition
1091 c write (iout,*) "After read_dist_constr nhpb",nhpb
1093 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1094 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1095 & modecalc.ne.10) then
1096 C If input structure hasn't been supplied from the PDB file read or generate
1098 if (iranconf.eq.0 .and. .not. extconf) then
1099 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1100 & write (iout,'(a)') 'Initial geometry will be read in.'
1102 read(inp,'(8f10.5)',end=36,err=36)
1103 & ((c(l,k),l=1,3),k=1,nres),
1104 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1105 write (iout,*) "Exit READ_CART"
1106 write (iout,'(8f10.5)')
1107 & ((c(l,k),l=1,3),k=1,nres),
1108 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1109 call int_from_cart1(.true.)
1110 write (iout,*) "Finish INT_TO_CART"
1113 dc(j,i)=c(j,i+1)-c(j,i)
1114 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1118 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1120 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1121 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1127 write (iout,*) "Calling read_ang"
1128 call read_angles(inp,*36)
1129 write (iout,*) "Calling chainbuild"
1130 call chainbuild_extconf
1133 36 write (iout,'(a)') 'Error reading angle file.'
1135 call mpi_finalize( MPI_COMM_WORLD,IERR )
1137 stop 'Error reading angle file.'
1139 else if (extconf) then
1140 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1141 & write (iout,'(a)') 'Extended chain initial geometry.'
1143 theta(i)=90d0*deg2rad
1146 phi(i)=180d0*deg2rad
1149 alph(i)=110d0*deg2rad
1152 omeg(i)=-120d0*deg2rad
1153 if (itype(i).le.0) omeg(i)=-omeg(i)
1155 c from old chainbuild
1157 C Define the origin and orientation of the coordinate system and locate the
1158 C first three CA's and SC(2).
1162 * Build the alpha-carbon chain.
1165 call locate_next_res(i)
1168 C First and last SC must coincide with the corresponding CA.
1172 dc_norm(j,nres+1)=0.0D0
1173 dc(j,nres+nres)=0.0D0
1174 dc_norm(j,nres+nres)=0.0D0
1176 c(j,nres+nres)=c(j,nres)
1179 C Define the origin and orientation of the coordinate system and locate the
1180 C first three CA's and SC(2).
1184 * Build the alpha-carbon chain.
1187 call locate_next_res(i)
1190 C First and last SC must coincide with the corresponding CA.
1194 dc_norm(j,nres+1)=0.0D0
1195 dc(j,nres+nres)=0.0D0
1196 dc_norm(j,nres+nres)=0.0D0
1198 c(j,nres+nres)=c(j,nres)
1203 if(me.eq.king.or..not.out1file)
1204 & write (iout,'(a)') 'Random-generated initial geometry.'
1208 if (me.eq.king .or. fg_rank.eq.0 .and. (
1209 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1213 call gen_rand_conf(itmp,*30)
1215 30 write (iout,*) 'Failed to generate random conformation',
1216 & ', itrial=',itrial
1217 write (*,*) 'Processor:',me,
1218 & ' Failed to generate random conformation',
1228 write (iout,'(a,i3,a)') 'Processor:',me,
1229 & ' error in generating random conformation.'
1230 write (*,'(a,i3,a)') 'Processor:',me,
1231 & ' error in generating random conformation.'
1234 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1239 & ' error in generating random conformation.'
1244 elseif (modecalc.eq.4) then
1245 read (inp,'(a)') intinname
1246 open (intin,file=intinname,status='old',err=333)
1247 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1248 & write (iout,'(a)') 'intinname',intinname
1249 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1251 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1253 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1255 stop 'Error opening angle file.'
1259 C Generate distance constraints, if the PDB structure is to be regularized.
1260 if (nthread.gt.0) then
1261 call read_threadbase
1264 if (me.eq.king .or. .not. out1file)
1266 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1267 write (iout,'(/a,i3,a)')
1268 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1269 write (iout,'(20i4)') (iss(i),i=1,ns)
1271 write(iout,*)"Running with dynamic disulfide-bond formation"
1273 write (iout,'(/a/)') 'Pre-formed links are:'
1279 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1280 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1286 if (ns.gt.0.and.dyn_ss) then
1290 forcon(i-nss)=forcon(i)
1297 dyn_ss_mask(iss(i))=.true.
1300 if (i2ndstr.gt.0) call secstrp2dihc
1301 c call geom_to_var(nvar,x)
1302 c call etotal(energia(0))
1303 c call enerprint(energia(0))
1304 c call briefout(0,etot)
1306 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1307 cd write (iout,'(a)') 'Variable list:'
1308 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1310 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1311 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1312 & 'Processor',myrank,': end reading molecular data.'
1316 c--------------------------------------------------------------------------
1317 logical function seq_comp(itypea,itypeb,length)
1319 integer length,itypea(length),itypeb(length)
1322 if (itypea(i).ne.itypeb(i)) then
1330 c-----------------------------------------------------------------------------
1331 subroutine read_bridge
1332 C Read information about disulfide bridges.
1333 implicit real*8 (a-h,o-z)
1334 include 'DIMENSIONS'
1338 include 'COMMON.IOUNITS'
1339 include 'COMMON.GEO'
1340 include 'COMMON.VAR'
1341 include 'COMMON.INTERACT'
1342 include 'COMMON.LOCAL'
1343 include 'COMMON.NAMES'
1344 include 'COMMON.CHAIN'
1345 include 'COMMON.FFIELD'
1346 include 'COMMON.SBRIDGE'
1347 include 'COMMON.HEADER'
1348 include 'COMMON.CONTROL'
1349 include 'COMMON.DBASE'
1350 include 'COMMON.THREAD'
1351 include 'COMMON.TIME1'
1352 include 'COMMON.SETUP'
1353 C Read bridging residues.
1354 read (inp,*) ns,(iss(i),i=1,ns)
1356 if(me.eq.king.or..not.out1file)
1357 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1358 C Check whether the specified bridging residues are cystines.
1360 if (itype(iss(i)).ne.1) then
1361 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1362 & 'Do you REALLY think that the residue ',
1363 & restyp(itype(iss(i))),i,
1364 & ' can form a disulfide bridge?!!!'
1365 write (*,'(2a,i3,a)')
1366 & 'Do you REALLY think that the residue ',
1367 & restyp(itype(iss(i))),i,
1368 & ' can form a disulfide bridge?!!!'
1370 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1375 C Read preformed bridges.
1377 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1379 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1382 C Check if the residues involved in bridges are in the specified list of
1383 C bridging residues.
1386 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1387 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1388 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1389 & ' contains residues present in other pairs.'
1390 write (*,'(a,i3,a)') 'Disulfide pair',i,
1391 & ' contains residues present in other pairs.'
1393 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1399 if (ihpb(i).eq.iss(j)) goto 10
1401 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1404 if (jhpb(i).eq.iss(j)) goto 20
1406 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1412 ihpb(i)=ihpb(i)+nres
1413 jhpb(i)=jhpb(i)+nres
1419 c----------------------------------------------------------------------------
1420 subroutine read_x(kanal,*)
1421 implicit real*8 (a-h,o-z)
1422 include 'DIMENSIONS'
1423 include 'COMMON.GEO'
1424 include 'COMMON.VAR'
1425 include 'COMMON.CHAIN'
1426 include 'COMMON.IOUNITS'
1427 include 'COMMON.CONTROL'
1428 include 'COMMON.LOCAL'
1429 include 'COMMON.INTERACT'
1430 c Read coordinates from input
1432 read(kanal,'(8f10.5)',end=10,err=10)
1433 & ((c(l,k),l=1,3),k=1,nres),
1434 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1437 c(j,2*nres)=c(j,nres)
1439 call int_from_cart1(.false.)
1442 dc(j,i)=c(j,i+1)-c(j,i)
1443 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1447 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1449 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1450 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1458 c----------------------------------------------------------------------------
1459 subroutine read_threadbase
1460 implicit real*8 (a-h,o-z)
1461 include 'DIMENSIONS'
1462 include 'COMMON.IOUNITS'
1463 include 'COMMON.GEO'
1464 include 'COMMON.VAR'
1465 include 'COMMON.INTERACT'
1466 include 'COMMON.LOCAL'
1467 include 'COMMON.NAMES'
1468 include 'COMMON.CHAIN'
1469 include 'COMMON.FFIELD'
1470 include 'COMMON.SBRIDGE'
1471 include 'COMMON.HEADER'
1472 include 'COMMON.CONTROL'
1473 include 'COMMON.DBASE'
1474 include 'COMMON.THREAD'
1475 include 'COMMON.TIME1'
1476 C Read pattern database for threading.
1477 read (icbase,*) nseq
1479 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1480 & nres_base(2,i),nres_base(3,i)
1481 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1483 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1484 c & nres_base(2,i),nres_base(3,i)
1485 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1489 if (weidis.eq.0.0D0) weidis=0.1D0
1498 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1499 write (iout,'(a,i5)') 'nexcl: ',nexcl
1500 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1503 c------------------------------------------------------------------------------
1504 subroutine setup_var
1505 implicit real*8 (a-h,o-z)
1506 include 'DIMENSIONS'
1507 include 'COMMON.IOUNITS'
1508 include 'COMMON.GEO'
1509 include 'COMMON.VAR'
1510 include 'COMMON.INTERACT'
1511 include 'COMMON.LOCAL'
1512 include 'COMMON.NAMES'
1513 include 'COMMON.CHAIN'
1514 include 'COMMON.FFIELD'
1515 include 'COMMON.SBRIDGE'
1516 include 'COMMON.HEADER'
1517 include 'COMMON.CONTROL'
1518 include 'COMMON.DBASE'
1519 include 'COMMON.THREAD'
1520 include 'COMMON.TIME1'
1521 C Set up variable list.
1527 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1529 ialph(i,1)=nvar+nside
1533 if (indphi.gt.0) then
1535 else if (indback.gt.0) then
1540 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1543 c----------------------------------------------------------------------------
1544 subroutine gen_dist_constr
1545 C Generate CA distance constraints.
1546 implicit real*8 (a-h,o-z)
1547 include 'DIMENSIONS'
1548 include 'COMMON.IOUNITS'
1549 include 'COMMON.GEO'
1550 include 'COMMON.VAR'
1551 include 'COMMON.INTERACT'
1552 include 'COMMON.LOCAL'
1553 include 'COMMON.NAMES'
1554 include 'COMMON.CHAIN'
1555 include 'COMMON.FFIELD'
1556 include 'COMMON.SBRIDGE'
1557 include 'COMMON.HEADER'
1558 include 'COMMON.CONTROL'
1559 include 'COMMON.DBASE'
1560 include 'COMMON.THREAD'
1561 include 'COMMON.TIME1'
1562 dimension itype_pdb(maxres)
1563 common /pizda/ itype_pdb
1565 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1566 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1567 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1569 do i=nstart_sup,nstart_sup+nsup-1
1570 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1571 cd & ' seq_pdb', restyp(itype_pdb(i))
1572 do j=i+2,nstart_sup+nsup-1
1574 ihpb(nhpb)=i+nstart_seq-nstart_sup
1575 jhpb(nhpb)=j+nstart_seq-nstart_sup
1577 dhpb(nhpb)=dist(i,j)
1580 cd write (iout,'(a)') 'Distance constraints:'
1585 cd if (ii.gt.nres) then
1590 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1591 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1592 cd & dhpb(i),forcon(i)
1596 c----------------------------------------------------------------------------
1598 implicit real*8 (a-h,o-z)
1599 include 'DIMENSIONS'
1600 include 'COMMON.MAP'
1601 include 'COMMON.IOUNITS'
1602 character*3 angid(4) /'THE','PHI','ALP','OME'/
1603 character*80 mapcard,ucase
1605 read (inp,'(a)') mapcard
1606 mapcard=ucase(mapcard)
1607 if (index(mapcard,'PHI').gt.0) then
1609 else if (index(mapcard,'THE').gt.0) then
1611 else if (index(mapcard,'ALP').gt.0) then
1613 else if (index(mapcard,'OME').gt.0) then
1616 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1617 stop 'Error - illegal variable spec in MAP card.'
1619 call readi (mapcard,'RES1',res1(imap),0)
1620 call readi (mapcard,'RES2',res2(imap),0)
1621 if (res1(imap).eq.0) then
1622 res1(imap)=res2(imap)
1623 else if (res2(imap).eq.0) then
1624 res2(imap)=res1(imap)
1626 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1628 & 'Error - illegal definition of variable group in MAP.'
1629 stop 'Error - illegal definition of variable group in MAP.'
1631 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1632 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1633 call readi(mapcard,'NSTEP',nstep(imap),0)
1634 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1636 & 'Illegal boundary and/or step size specification in MAP.'
1637 stop 'Illegal boundary and/or step size specification in MAP.'
1642 c----------------------------------------------------------------------------
1644 implicit real*8 (a-h,o-z)
1645 include 'DIMENSIONS'
1646 include 'COMMON.IOUNITS'
1647 include 'COMMON.GEO'
1648 include 'COMMON.CSA'
1649 include 'COMMON.BANK'
1650 include 'COMMON.CONTROL'
1652 character*620 mcmcard
1653 call card_concat(mcmcard)
1655 call readi(mcmcard,'NCONF',nconf,50)
1656 call readi(mcmcard,'NADD',nadd,0)
1657 call readi(mcmcard,'JSTART',jstart,1)
1658 call readi(mcmcard,'JEND',jend,1)
1659 call readi(mcmcard,'NSTMAX',nstmax,500000)
1660 call readi(mcmcard,'N0',n0,1)
1661 call readi(mcmcard,'N1',n1,6)
1662 call readi(mcmcard,'N2',n2,4)
1663 call readi(mcmcard,'N3',n3,0)
1664 call readi(mcmcard,'N4',n4,0)
1665 call readi(mcmcard,'N5',n5,0)
1666 call readi(mcmcard,'N6',n6,10)
1667 call readi(mcmcard,'N7',n7,0)
1668 call readi(mcmcard,'N8',n8,0)
1669 call readi(mcmcard,'N9',n9,0)
1670 call readi(mcmcard,'N14',n14,0)
1671 call readi(mcmcard,'N15',n15,0)
1672 call readi(mcmcard,'N16',n16,0)
1673 call readi(mcmcard,'N17',n17,0)
1674 call readi(mcmcard,'N18',n18,0)
1676 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1678 call readi(mcmcard,'NDIFF',ndiff,2)
1679 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1680 call readi(mcmcard,'IS1',is1,1)
1681 call readi(mcmcard,'IS2',is2,8)
1682 call readi(mcmcard,'NRAN0',nran0,4)
1683 call readi(mcmcard,'NRAN1',nran1,2)
1684 call readi(mcmcard,'IRR',irr,1)
1685 call readi(mcmcard,'NSEED',nseed,20)
1686 call readi(mcmcard,'NTOTAL',ntotal,10000)
1687 call reada(mcmcard,'CUT1',cut1,2.0d0)
1688 call reada(mcmcard,'CUT2',cut2,5.0d0)
1689 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1690 call readi(mcmcard,'ICMAX',icmax,3)
1691 call readi(mcmcard,'IRESTART',irestart,0)
1692 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1695 call reada(mcmcard,'DELE',dele,20.0d0)
1696 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1697 call readi(mcmcard,'IREF',iref,0)
1698 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1699 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1700 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1701 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1702 write (iout,*) "NCONF_IN",nconf_in
1705 c----------------------------------------------------------------------------
1706 cfmc subroutine mcmfread
1707 cfmc implicit real*8 (a-h,o-z)
1708 cfmc include 'DIMENSIONS'
1709 cfmc include 'COMMON.MCMF'
1710 cfmc include 'COMMON.IOUNITS'
1711 cfmc include 'COMMON.GEO'
1712 cfmc character*80 ucase
1713 cfmc character*620 mcmcard
1714 cfmc call card_concat(mcmcard)
1716 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1717 cfmc write(iout,*)'MAXRANT=',maxrant
1718 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1719 cfmc write(iout,*)'MAXFAM=',maxfam
1720 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1721 cfmc write(iout,*)'NNET1=',nnet1
1722 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1723 cfmc write(iout,*)'NNET2=',nnet2
1724 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1725 cfmc write(iout,*)'NNET3=',nnet3
1726 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1727 cfmc write(iout,*)'ILASTT=',ilastt
1728 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1729 cfmc write(iout,*)'MAXSTR=',maxstr
1730 cfmc maxstr_f=maxstr/maxfam
1731 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1732 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1733 cfmc write(iout,*)'NMCMF=',nmcmf
1734 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1735 cfmc write(iout,*)'IFOCUS=',ifocus
1736 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1737 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1738 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1739 cfmc write(iout,*)'INTPRT=',intprt
1740 cfmc call readi(mcmcard,'IPRT',iprt,100)
1741 cfmc write(iout,*)'IPRT=',iprt
1742 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1743 cfmc write(iout,*)'IMAXTR=',imaxtr
1744 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1745 cfmc write(iout,*)'MAXEVEN=',maxeven
1746 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1747 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1748 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1749 cfmc write(iout,*)'INIMIN=',inimin
1750 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1751 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1752 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1753 cfmc write(iout,*)'NTHREAD=',nthread
1754 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1755 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1756 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1757 cfmc write(iout,*)'MAXPERT=',maxpert
1758 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1759 cfmc write(iout,*)'IRMSD=',irmsd
1760 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1761 cfmc write(iout,*)'DENEMIN=',denemin
1762 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1763 cfmc write(iout,*)'RCUT1S=',rcut1s
1764 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1765 cfmc write(iout,*)'RCUT1E=',rcut1e
1766 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1767 cfmc write(iout,*)'RCUT2S=',rcut2s
1768 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1769 cfmc write(iout,*)'RCUT2E=',rcut2e
1770 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1771 cfmc write(iout,*)'DPERT1=',d_pert1
1772 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1773 cfmc write(iout,*)'DPERT1A=',d_pert1a
1774 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1775 cfmc write(iout,*)'DPERT2=',d_pert2
1776 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1777 cfmc write(iout,*)'DPERT2A=',d_pert2a
1778 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1779 cfmc write(iout,*)'DPERT2B=',d_pert2b
1780 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1781 cfmc write(iout,*)'DPERT2C=',d_pert2c
1782 cfmc d_pert1=deg2rad*d_pert1
1783 cfmc d_pert1a=deg2rad*d_pert1a
1784 cfmc d_pert2=deg2rad*d_pert2
1785 cfmc d_pert2a=deg2rad*d_pert2a
1786 cfmc d_pert2b=deg2rad*d_pert2b
1787 cfmc d_pert2c=deg2rad*d_pert2c
1788 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1789 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1790 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1791 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1792 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1793 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1794 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1795 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1796 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1797 cfmc write(iout,*)'RCUTINI=',rcutini
1798 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1799 cfmc write(iout,*)'GRAT=',grat
1800 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1801 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1805 c----------------------------------------------------------------------------
1807 implicit real*8 (a-h,o-z)
1808 include 'DIMENSIONS'
1809 include 'COMMON.MCM'
1810 include 'COMMON.MCE'
1811 include 'COMMON.IOUNITS'
1813 character*320 mcmcard
1814 call card_concat(mcmcard)
1815 call readi(mcmcard,'MAXACC',maxacc,100)
1816 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1817 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1818 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1819 call readi(mcmcard,'MAXREPM',maxrepm,200)
1820 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1821 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1822 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1823 call reada(mcmcard,'E_UP',e_up,5.0D0)
1824 call reada(mcmcard,'DELTE',delte,0.1D0)
1825 call readi(mcmcard,'NSWEEP',nsweep,5)
1826 call readi(mcmcard,'NSTEPH',nsteph,0)
1827 call readi(mcmcard,'NSTEPC',nstepc,0)
1828 call reada(mcmcard,'TMIN',tmin,298.0D0)
1829 call reada(mcmcard,'TMAX',tmax,298.0D0)
1830 call readi(mcmcard,'NWINDOW',nwindow,0)
1831 call readi(mcmcard,'PRINT_MC',print_mc,0)
1832 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1833 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1834 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1835 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1836 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1837 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1838 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1839 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1840 if (nwindow.gt.0) then
1841 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1843 winlen(i)=winend(i)-winstart(i)+1
1846 if (tmax.lt.tmin) tmax=tmin
1847 if (tmax.eq.tmin) then
1851 if (nstepc.gt.0 .and. nsteph.gt.0) then
1852 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1853 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1855 C Probabilities of different move types
1856 sumpro_type(0)=0.0D0
1857 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1858 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1859 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1860 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1861 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1862 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1863 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1865 print *,'i',i,' sumprotype',sumpro_type(i)
1866 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1867 print *,'i',i,' sumprotype',sumpro_type(i)
1871 c----------------------------------------------------------------------------
1872 subroutine read_minim
1873 implicit real*8 (a-h,o-z)
1874 include 'DIMENSIONS'
1875 include 'COMMON.MINIM'
1876 include 'COMMON.IOUNITS'
1877 include 'COMMON.CONTROL'
1878 include 'COMMON.SETUP'
1880 character*320 minimcard
1881 call card_concat(minimcard)
1882 call readi(minimcard,'MAXMIN',maxmin,2000)
1883 call readi(minimcard,'MAXFUN',maxfun,5000)
1884 call readi(minimcard,'MINMIN',minmin,maxmin)
1885 call readi(minimcard,'MINFUN',minfun,maxmin)
1886 call reada(minimcard,'TOLF',tolf,1.0D-2)
1887 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1888 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1889 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1890 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1892 if (.not. out1file .or. me.eq.king) then
1894 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1895 & 'Options in energy minimization:'
1896 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1897 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1898 & 'MinMin:',MinMin,' MinFun:',MinFun,
1899 & ' TolF:',TolF,' RTolF:',RTolF
1905 c----------------------------------------------------------------------------
1906 subroutine read_angles(kanal,*)
1907 implicit real*8 (a-h,o-z)
1908 include 'DIMENSIONS'
1909 include 'COMMON.GEO'
1910 include 'COMMON.VAR'
1911 include 'COMMON.CHAIN'
1912 include 'COMMON.IOUNITS'
1913 include 'COMMON.CONTROL'
1914 c Read angles from input
1916 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1917 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1918 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1919 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1922 c 9/7/01 avoid 180 deg valence angle
1923 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1925 theta(i)=deg2rad*theta(i)
1926 phi(i)=deg2rad*phi(i)
1927 alph(i)=deg2rad*alph(i)
1928 omeg(i)=deg2rad*omeg(i)
1933 c----------------------------------------------------------------------------
1934 subroutine reada(rekord,lancuch,wartosc,default)
1936 character*(*) rekord,lancuch
1937 double precision wartosc,default
1940 iread=index(rekord,lancuch)
1941 if (iread.eq.0) then
1945 iread=iread+ilen(lancuch)+1
1946 read (rekord(iread:),*,err=10,end=10) wartosc
1951 c----------------------------------------------------------------------------
1952 subroutine readi(rekord,lancuch,wartosc,default)
1954 character*(*) rekord,lancuch
1955 integer wartosc,default
1958 iread=index(rekord,lancuch)
1959 if (iread.eq.0) then
1963 iread=iread+ilen(lancuch)+1
1964 read (rekord(iread:),*,err=10,end=10) wartosc
1969 c----------------------------------------------------------------------------
1970 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1973 integer tablica(dim),default
1974 character*(*) rekord,lancuch
1981 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1982 if (iread.eq.0) return
1983 iread=iread+ilen(lancuch)+1
1984 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1987 c----------------------------------------------------------------------------
1988 subroutine multreada(rekord,lancuch,tablica,dim,default)
1991 double precision tablica(dim),default
1992 character*(*) rekord,lancuch
1999 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2000 if (iread.eq.0) return
2001 iread=iread+ilen(lancuch)+1
2002 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2005 c----------------------------------------------------------------------------
2006 subroutine openunits
2007 implicit real*8 (a-h,o-z)
2008 include 'DIMENSIONS'
2011 character*16 form,nodename
2014 include 'COMMON.SETUP'
2015 include 'COMMON.IOUNITS'
2017 include 'COMMON.CONTROL'
2018 integer lenpre,lenpot,ilen,lentmp
2020 character*3 out1file_text,ucase
2023 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2024 call getenv_loc("PREFIX",prefix)
2026 call getenv_loc("POT",pot)
2027 call getenv_loc("DIRTMP",tmpdir)
2028 call getenv_loc("CURDIR",curdir)
2029 call getenv_loc("OUT1FILE",out1file_text)
2030 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2031 out1file_text=ucase(out1file_text)
2032 if (out1file_text(1:1).eq."Y") then
2035 out1file=fg_rank.gt.0
2040 if (lentmp.gt.0) then
2041 write (*,'(80(1h!))')
2042 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2043 write (*,'(80(1h!))')
2044 write (*,*)"All output files will be on node /tmp directory."
2046 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2047 if (me.eq.king) then
2048 write (*,*) "The master node is ",nodename
2049 else if (fg_rank.eq.0) then
2050 write (*,*) "I am the CG slave node ",nodename
2052 write (*,*) "I am the FG slave node ",nodename
2055 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2056 lenpre = lentmp+lenpre+1
2058 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2059 C Get the names and open the input files
2060 #if defined(WINIFL) || defined(WINPGI)
2061 open(1,file=pref_orig(:ilen(pref_orig))//
2062 & '.inp',status='old',readonly,shared)
2063 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2064 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2065 C Get parameter filenames and open the parameter files.
2066 call getenv_loc('BONDPAR',bondname)
2067 open (ibond,file=bondname,status='old',readonly,shared)
2068 call getenv_loc('THETPAR',thetname)
2069 open (ithep,file=thetname,status='old',readonly,shared)
2071 call getenv_loc('THETPARPDB',thetname_pdb)
2072 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2074 call getenv_loc('ROTPAR',rotname)
2075 open (irotam,file=rotname,status='old',readonly,shared)
2077 call getenv_loc('ROTPARPDB',rotname_pdb)
2078 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2080 call getenv_loc('TORPAR',torname)
2081 open (itorp,file=torname,status='old',readonly,shared)
2082 call getenv_loc('TORDPAR',tordname)
2083 open (itordp,file=tordname,status='old',readonly,shared)
2084 call getenv_loc('FOURIER',fouriername)
2085 open (ifourier,file=fouriername,status='old',readonly,shared)
2086 call getenv_loc('ELEPAR',elename)
2087 open (ielep,file=elename,status='old',readonly,shared)
2088 call getenv_loc('SIDEPAR',sidename)
2089 open (isidep,file=sidename,status='old',readonly,shared)
2090 #elif (defined CRAY) || (defined AIX)
2091 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2093 c print *,"Processor",myrank," opened file 1"
2094 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2095 c print *,"Processor",myrank," opened file 9"
2096 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2097 C Get parameter filenames and open the parameter files.
2098 call getenv_loc('BONDPAR',bondname)
2099 open (ibond,file=bondname,status='old',action='read')
2100 c print *,"Processor",myrank," opened file IBOND"
2101 call getenv_loc('THETPAR',thetname)
2102 open (ithep,file=thetname,status='old',action='read')
2103 c print *,"Processor",myrank," opened file ITHEP"
2105 call getenv_loc('THETPARPDB',thetname_pdb)
2106 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2108 call getenv_loc('ROTPAR',rotname)
2109 open (irotam,file=rotname,status='old',action='read')
2110 c print *,"Processor",myrank," opened file IROTAM"
2112 call getenv_loc('ROTPARPDB',rotname_pdb)
2113 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2115 call getenv_loc('TORPAR',torname)
2116 open (itorp,file=torname,status='old',action='read')
2117 c print *,"Processor",myrank," opened file ITORP"
2118 call getenv_loc('TORDPAR',tordname)
2119 open (itordp,file=tordname,status='old',action='read')
2120 c print *,"Processor",myrank," opened file ITORDP"
2121 call getenv_loc('SCCORPAR',sccorname)
2122 open (isccor,file=sccorname,status='old',action='read')
2123 c print *,"Processor",myrank," opened file ISCCOR"
2124 call getenv_loc('FOURIER',fouriername)
2125 open (ifourier,file=fouriername,status='old',action='read')
2126 c print *,"Processor",myrank," opened file IFOURIER"
2127 call getenv_loc('ELEPAR',elename)
2128 open (ielep,file=elename,status='old',action='read')
2129 c print *,"Processor",myrank," opened file IELEP"
2130 call getenv_loc('SIDEPAR',sidename)
2131 open (isidep,file=sidename,status='old',action='read')
2132 c print *,"Processor",myrank," opened file ISIDEP"
2133 c print *,"Processor",myrank," opened parameter files"
2135 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2136 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2137 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2138 C Get parameter filenames and open the parameter files.
2139 call getenv_loc('BONDPAR',bondname)
2140 open (ibond,file=bondname,status='old')
2141 call getenv_loc('THETPAR',thetname)
2142 open (ithep,file=thetname,status='old')
2144 call getenv_loc('THETPARPDB',thetname_pdb)
2145 open (ithep_pdb,file=thetname_pdb,status='old')
2147 call getenv_loc('ROTPAR',rotname)
2148 open (irotam,file=rotname,status='old')
2150 call getenv_loc('ROTPARPDB',rotname_pdb)
2151 open (irotam_pdb,file=rotname_pdb,status='old')
2153 call getenv_loc('TORPAR',torname)
2154 open (itorp,file=torname,status='old')
2155 call getenv_loc('TORDPAR',tordname)
2156 open (itordp,file=tordname,status='old')
2157 call getenv_loc('SCCORPAR',sccorname)
2158 open (isccor,file=sccorname,status='old')
2159 call getenv_loc('FOURIER',fouriername)
2160 open (ifourier,file=fouriername,status='old')
2161 call getenv_loc('ELEPAR',elename)
2162 open (ielep,file=elename,status='old')
2163 call getenv_loc('SIDEPAR',sidename)
2164 open (isidep,file=sidename,status='old')
2165 call getenv_loc('LIPTRANPAR',liptranname)
2166 open (iliptranpar,file=liptranname,status='old')
2168 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2170 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2171 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2172 C Get parameter filenames and open the parameter files.
2173 call getenv_loc('BONDPAR',bondname)
2174 open (ibond,file=bondname,status='old',readonly)
2175 call getenv_loc('THETPAR',thetname)
2176 open (ithep,file=thetname,status='old',readonly)
2177 call getenv_loc('ROTPAR',rotname)
2178 open (irotam,file=rotname,status='old',readonly)
2179 call getenv_loc('TORPAR',torname)
2180 open (itorp,file=torname,status='old',readonly)
2181 call getenv_loc('TORDPAR',tordname)
2182 open (itordp,file=tordname,status='old',readonly)
2183 call getenv_loc('SCCORPAR',sccorname)
2184 open (isccor,file=sccorname,status='old',readonly)
2186 call getenv_loc('THETPARPDB',thetname_pdb)
2187 c print *,"thetname_pdb ",thetname_pdb
2188 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2189 c print *,ithep_pdb," opened"
2191 call getenv_loc('FOURIER',fouriername)
2192 open (ifourier,file=fouriername,status='old',readonly)
2193 call getenv_loc('ELEPAR',elename)
2194 open (ielep,file=elename,status='old',readonly)
2195 call getenv_loc('SIDEPAR',sidename)
2196 open (isidep,file=sidename,status='old',readonly)
2197 call getenv_loc('LIPTRANPAR',liptranname)
2198 open (iliptranpar,file=liptranname,status='old',readonly)
2200 call getenv_loc('ROTPARPDB',rotname_pdb)
2201 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2206 C 8/9/01 In the newest version SCp interaction constants are read from a file
2207 C Use -DOLDSCP to use hard-coded constants instead.
2209 call getenv_loc('SCPPAR',scpname)
2210 #if defined(WINIFL) || defined(WINPGI)
2211 open (iscpp,file=scpname,status='old',readonly,shared)
2212 #elif (defined CRAY) || (defined AIX)
2213 open (iscpp,file=scpname,status='old',action='read')
2215 open (iscpp,file=scpname,status='old')
2217 open (iscpp,file=scpname,status='old',readonly)
2220 call getenv_loc('PATTERN',patname)
2221 #if defined(WINIFL) || defined(WINPGI)
2222 open (icbase,file=patname,status='old',readonly,shared)
2223 #elif (defined CRAY) || (defined AIX)
2224 open (icbase,file=patname,status='old',action='read')
2226 open (icbase,file=patname,status='old')
2228 open (icbase,file=patname,status='old',readonly)
2231 C Open output file only for CG processes
2232 c print *,"Processor",myrank," fg_rank",fg_rank
2233 if (fg_rank.eq.0) then
2235 if (nodes.eq.1) then
2238 npos = dlog10(dfloat(nodes-1))+1
2240 if (npos.lt.3) npos=3
2241 write (liczba,'(i1)') npos
2242 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2244 write (liczba,form) me
2245 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2246 & liczba(:ilen(liczba))
2247 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2249 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2251 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2252 & liczba(:ilen(liczba))//'.mol2'
2253 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2254 & liczba(:ilen(liczba))//'.stat'
2256 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2257 & //liczba(:ilen(liczba))//'.stat')
2258 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2261 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2262 & liczba(:ilen(liczba))//'.const'
2267 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2268 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2269 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2270 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2271 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2273 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2275 rest2name=prefix(:ilen(prefix))//'.rst'
2277 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2280 #if defined(AIX) || defined(PGI)
2281 if (me.eq.king .or. .not. out1file)
2282 & open(iout,file=outname,status='unknown')
2284 if (fg_rank.gt.0) then
2285 write (liczba,'(i3.3)') myrank/nfgtasks
2286 write (ll,'(bz,i3.3)') fg_rank
2287 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2292 open(igeom,file=intname,status='unknown',position='append')
2293 open(ipdb,file=pdbname,status='unknown')
2294 open(imol2,file=mol2name,status='unknown')
2295 open(istat,file=statname,status='unknown',position='append')
2297 c1out open(iout,file=outname,status='unknown')
2300 if (me.eq.king .or. .not.out1file)
2301 & open(iout,file=outname,status='unknown')
2303 if (fg_rank.gt.0) then
2304 write (liczba,'(i3.3)') myrank/nfgtasks
2305 write (ll,'(bz,i3.3)') fg_rank
2306 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2311 open(igeom,file=intname,status='unknown',access='append')
2312 open(ipdb,file=pdbname,status='unknown')
2313 open(imol2,file=mol2name,status='unknown')
2314 open(istat,file=statname,status='unknown',access='append')
2316 c1out open(iout,file=outname,status='unknown')
2319 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2320 csa_seed=prefix(:lenpre)//'.CSA.seed'
2321 csa_history=prefix(:lenpre)//'.CSA.history'
2322 csa_bank=prefix(:lenpre)//'.CSA.bank'
2323 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2324 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2325 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2326 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2327 csa_int=prefix(:lenpre)//'.int'
2328 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2329 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2330 csa_in=prefix(:lenpre)//'.CSA.in'
2331 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2334 write (iout,'(80(1h-))')
2335 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2336 write (iout,'(80(1h-))')
2337 write (iout,*) "Input file : ",
2338 & pref_orig(:ilen(pref_orig))//'.inp'
2339 write (iout,*) "Output file : ",
2340 & outname(:ilen(outname))
2342 write (iout,*) "Sidechain potential file : ",
2343 & sidename(:ilen(sidename))
2345 write (iout,*) "SCp potential file : ",
2346 & scpname(:ilen(scpname))
2348 write (iout,*) "Electrostatic potential file : ",
2349 & elename(:ilen(elename))
2350 write (iout,*) "Cumulant coefficient file : ",
2351 & fouriername(:ilen(fouriername))
2352 write (iout,*) "Torsional parameter file : ",
2353 & torname(:ilen(torname))
2354 write (iout,*) "Double torsional parameter file : ",
2355 & tordname(:ilen(tordname))
2356 write (iout,*) "SCCOR parameter file : ",
2357 & sccorname(:ilen(sccorname))
2358 write (iout,*) "Bond & inertia constant file : ",
2359 & bondname(:ilen(bondname))
2360 write (iout,*) "Bending parameter file : ",
2361 & thetname(:ilen(thetname))
2362 write (iout,*) "Rotamer parameter file : ",
2363 & rotname(:ilen(rotname))
2364 write (iout,*) "Thetpdb parameter file : ",
2365 & thetname_pdb(:ilen(thetname_pdb))
2366 write (iout,*) "Threading database : ",
2367 & patname(:ilen(patname))
2369 &write (iout,*)" DIRTMP : ",
2371 write (iout,'(80(1h-))')
2375 c----------------------------------------------------------------------------
2376 subroutine card_concat(card)
2377 implicit real*8 (a-h,o-z)
2378 include 'DIMENSIONS'
2379 include 'COMMON.IOUNITS'
2381 character*80 karta,ucase
2383 read (inp,'(a)') karta
2386 do while (karta(80:80).eq.'&')
2387 card=card(:ilen(card)+1)//karta(:79)
2388 read (inp,'(a)') karta
2391 card=card(:ilen(card)+1)//karta
2394 c----------------------------------------------------------------------------------
2396 implicit real*8 (a-h,o-z)
2397 include 'DIMENSIONS'
2398 include 'COMMON.CHAIN'
2399 include 'COMMON.IOUNITS'
2401 include 'COMMON.CONTROL'
2402 open(irest2,file=rest2name,status='unknown')
2403 read(irest2,*) totT,EK,potE,totE,t_bath
2406 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2409 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2411 if(usampl.or.homol_nset.gt.1) then
2412 read (irest2,*) iset
2417 c---------------------------------------------------------------------------------
2418 subroutine read_fragments
2419 implicit real*8 (a-h,o-z)
2420 include 'DIMENSIONS'
2424 include 'COMMON.SETUP'
2425 include 'COMMON.CHAIN'
2426 include 'COMMON.IOUNITS'
2428 include 'COMMON.CONTROL'
2430 read(inp,*) nset,nfrag,npair,nfrag_back
2431 if(me.eq.king.or..not.out1file)
2432 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2433 & " nfrag_back",nfrag_back
2435 read(inp,*) mset(iset1)
2437 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2439 if(me.eq.king.or..not.out1file)
2440 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2441 & ifrag(2,i,iset1), qinfrag(i,iset1)
2444 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2446 if(me.eq.king.or..not.out1file)
2447 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2448 & ipair(2,i,iset1), qinpair(i,iset1)
2451 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2452 & wfrag_back(3,i,iset1),
2453 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2454 if(me.eq.king.or..not.out1file)
2455 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2456 & wfrag_back(2,i,iset1),
2457 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2458 & ifrag_back(2,i,iset1)
2463 C---------------------------------------------------------------------------
2464 subroutine read_afminp
2465 implicit real*8 (a-h,o-z)
2466 include 'DIMENSIONS'
2470 include 'COMMON.SETUP'
2471 include 'COMMON.CONTROL'
2472 include 'COMMON.CHAIN'
2473 include 'COMMON.IOUNITS'
2474 include 'COMMON.SBRIDGE'
2475 character*320 afmcard
2476 c print *, "wchodze"
2477 call card_concat(afmcard)
2478 call readi(afmcard,"BEG",afmbeg,0)
2479 call readi(afmcard,"END",afmend,0)
2480 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2481 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2482 print *,'FORCE=' ,forceAFMconst
2483 CCCC NOW PROPERTIES FOR AFM
2486 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2488 distafminit=dsqrt(distafminit)
2489 print *,'initdist',distafminit
2493 c-------------------------------------------------------------------------------
2494 subroutine read_saxs_constr
2495 implicit real*8 (a-h,o-z)
2496 include 'DIMENSIONS'
2500 include 'COMMON.SETUP'
2501 include 'COMMON.CONTROL'
2502 include 'COMMON.CHAIN'
2503 include 'COMMON.IOUNITS'
2504 include 'COMMON.SBRIDGE'
2505 double precision cm(3)
2507 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2509 if (saxs_mode.eq.0) then
2510 c SAXS distance distribution
2512 read(inp,*) distsaxs(i),Psaxs(i)
2516 Cnorm = Cnorm + Psaxs(i)
2518 write (iout,*) "Cnorm",Cnorm
2520 Psaxs(i)=Psaxs(i)/Cnorm
2522 write (iout,*) "Normalized distance distribution from SAXS"
2524 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2528 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2530 write (iout,*) "Wsaxs0",Wsaxs0
2534 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2541 cm(j)=cm(j)+Csaxs(j,i)
2549 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2552 write (iout,*) "SAXS sphere coordinates"
2554 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2559 c-------------------------------------------------------------------------------
2560 subroutine read_dist_constr
2561 implicit real*8 (a-h,o-z)
2562 include 'DIMENSIONS'
2566 include 'COMMON.SETUP'
2567 include 'COMMON.CONTROL'
2568 include 'COMMON.CHAIN'
2569 include 'COMMON.IOUNITS'
2570 include 'COMMON.SBRIDGE'
2571 integer ifrag_(2,100),ipair_(2,1000)
2572 double precision wfrag_(100),wpair_(1000)
2573 character*5000 controlcard
2574 logical normalize,next
2576 double precision xlink(4,0:4) /
2578 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2579 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2580 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2581 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2582 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2583 c print *, "WCHODZE"
2584 write (iout,*) "Calling read_dist_constr"
2585 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2595 call card_concat(controlcard)
2596 next = index(controlcard,"NEXT").gt.0
2597 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2598 write (iout,*) "restr_type",restr_type
2599 call readi(controlcard,"NFRAG",nfrag_,0)
2600 call readi(controlcard,"NFRAG",nfrag_,0)
2601 call readi(controlcard,"NPAIR",npair_,0)
2602 call readi(controlcard,"NDIST",ndist_,0)
2603 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2604 if (restr_type.eq.10)
2605 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2606 if (restr_type.eq.12)
2607 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2608 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2609 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2610 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2611 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2612 normalize = index(controlcard,"NORMALIZE").gt.0
2613 write (iout,*) "WBOLTZD",wboltzd
2614 write (iout,*) "SCAL_PEAK",scal_peak
2615 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2616 write (iout,*) "IFRAG"
2618 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2620 write (iout,*) "IPAIR"
2622 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2624 if (nfrag_.gt.0) write (iout,*)
2625 & "Distance restraints as generated from reference structure"
2627 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2628 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2629 & ifrag_(2,i)=nstart_sup+nsup-1
2630 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2632 if (wfrag_(i).eq.0.0d0) cycle
2633 do j=ifrag_(1,i),ifrag_(2,i)-1
2634 do k=j+1,ifrag_(2,i)
2635 c write (iout,*) "j",j," k",k
2637 if (restr_type.eq.1) then
2643 forcon(nhpb)=wfrag_(i)
2644 else if (constr_dist.eq.2) then
2645 if (ddjk.le.dist_cut) then
2651 forcon(nhpb)=wfrag_(i)
2653 else if (restr_type.eq.3) then
2659 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2662 if (.not.out1file .or. me.eq.king)
2663 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2664 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2666 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2667 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2673 if (wpair_(i).eq.0.0d0) cycle
2681 do j=ifrag_(1,ii),ifrag_(2,ii)
2682 do k=ifrag_(1,jj),ifrag_(2,jj)
2684 if (restr_type.eq.1) then
2690 forcon(nhpb)=wpair_(i)
2691 else if (constr_dist.eq.2) then
2692 if (ddjk.le.dist_cut) then
2698 forcon(nhpb)=wpair_(i)
2700 else if (restr_type.eq.3) then
2706 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2709 if (.not.out1file .or. me.eq.king)
2710 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2711 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2713 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2714 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2721 write (iout,*) "Distance restraints as read from input"
2723 if (restr_type.eq.12) then
2724 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2725 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2726 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2727 & fordepth_peak(nhpb_peak+1),npeak
2728 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2729 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2730 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2731 c & fordepth_peak(nhpb_peak+1),npeak
2732 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2733 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2734 nhpb_peak=nhpb_peak+1
2735 irestr_type_peak(nhpb_peak)=12
2736 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2739 if (.not.out1file .or. me.eq.king)
2740 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2741 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2742 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2743 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2744 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2746 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2747 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2748 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2749 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2750 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2752 if (ibecarb_peak(nhpb_peak).gt.0) then
2753 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2754 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2756 else if (restr_type.eq.11) then
2757 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2758 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2759 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2760 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2762 irestr_type(nhpb)=11
2764 if (.not.out1file .or. me.eq.king)
2765 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2766 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2767 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2769 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2770 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2771 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2773 if (ibecarb(nhpb).gt.0) then
2774 ihpb(nhpb)=ihpb(nhpb)+nres
2775 jhpb(nhpb)=jhpb(nhpb)+nres
2777 else if (restr_type.eq.10) then
2778 c Cross-lonk Markov-like potential
2779 call card_concat(controlcard)
2780 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2781 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2783 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2784 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2785 if (index(controlcard,"ZL").gt.0) then
2787 else if (index(controlcard,"ADH").gt.0) then
2789 else if (index(controlcard,"PDH").gt.0) then
2791 else if (index(controlcard,"DSS").gt.0) then
2796 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2797 & xlink(1,link_type))
2798 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2799 & xlink(2,link_type))
2800 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2801 & xlink(3,link_type))
2802 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2803 & xlink(4,link_type))
2804 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2805 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2806 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2807 if (forcon(nhpb+1).le.0.0d0 .or.
2808 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2810 irestr_type(nhpb)=10
2811 if (ibecarb(nhpb).gt.0) then
2812 ihpb(nhpb)=ihpb(nhpb)+nres
2813 jhpb(nhpb)=jhpb(nhpb)+nres
2816 if (.not.out1file .or. me.eq.king)
2817 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2818 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2819 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2822 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2823 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2824 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2829 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2830 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2831 if (forcon(nhpb+1).gt.0.0d0) then
2833 if (dhpb1(nhpb).eq.0.0d0) then
2838 if (ibecarb(nhpb).gt.0) then
2839 ihpb(nhpb)=ihpb(nhpb)+nres
2840 jhpb(nhpb)=jhpb(nhpb)+nres
2842 if (dhpb(nhpb).eq.0.0d0)
2843 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2846 if (.not.out1file .or. me.eq.king)
2847 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2848 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2850 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2851 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2854 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2855 C if (forcon(nhpb+1).gt.0.0d0) then
2857 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2865 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2866 & fordepthmax=fordepth(i)
2869 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2874 c-------------------------------------------------------------------------------
2876 subroutine read_constr_homology
2878 include 'DIMENSIONS'
2882 include 'COMMON.SETUP'
2883 include 'COMMON.CONTROL'
2884 include 'COMMON.CHAIN'
2885 include 'COMMON.IOUNITS'
2887 include 'COMMON.GEO'
2888 include 'COMMON.INTERACT'
2889 include 'COMMON.NAMES'
2891 c For new homol impl
2893 include 'COMMON.VAR'
2896 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2898 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2899 c & sigma_odl_temp(maxres,maxres,max_template)
2901 character*24 model_ki_dist, model_ki_angle
2902 character*500 controlcard
2903 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2904 logical lprn /.true./
2909 c FP - Nov. 2014 Temporary specifications for new vars
2911 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2913 double precision, dimension (max_template,maxres) :: rescore
2914 double precision, dimension (max_template,maxres) :: rescore2
2915 double precision, dimension (max_template,maxres) :: rescore3
2916 character*24 pdbfile,tpl_k_rescore
2917 c -----------------------------------------------------------------
2918 c Reading multiple PDB ref structures and calculation of retraints
2919 c not using pre-computed ones stored in files model_ki_{dist,angle}
2921 c -----------------------------------------------------------------
2924 c Alternative: reading from input
2925 call card_concat(controlcard)
2926 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2927 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2928 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2929 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2930 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2931 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2932 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2933 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2934 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2935 if(.not.read2sigma.and.start_from_model) then
2936 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2937 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2938 start_from_model=.false.
2940 if(start_from_model .and. (me.eq.king .or. .not. out1file))
2941 & write(iout,*) 'START_FROM_MODELS is ON'
2942 if(start_from_model .and. rest) then
2943 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2944 write(iout,*) 'START_FROM_MODELS is OFF'
2945 write(iout,*) 'remove restart keyword from input'
2948 if (homol_nset.gt.1)then
2949 call card_concat(controlcard)
2950 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2951 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2952 write(iout,*) "iset homology_weight "
2954 write(iout,*) i,waga_homology(i)
2957 iset=mod(kolor,homol_nset)+1
2960 waga_homology(1)=1.0
2963 cd write (iout,*) "nnt",nnt," nct",nct
2970 c write(iout,*) 'nnt=',nnt,'nct=',nct
2973 do k=1,constr_homology
2986 if (read_homol_frag) then
2987 call read_klapaucjusz
2990 do k=1,constr_homology
2992 read(inp,'(a)') pdbfile
2993 if(me.eq.king .or. .not. out1file)
2994 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2995 & pdbfile(:ilen(pdbfile))
2996 open(ipdbin,file=pdbfile,status='old',err=33)
2998 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2999 & pdbfile(:ilen(pdbfile))
3002 c print *,'Begin reading pdb data'
3004 c Files containing res sim or local scores (former containing sigmas)
3007 write(kic2,'(bz,i2.2)') k
3009 tpl_k_rescore="template"//kic2//".sco"
3012 if (read2sigma) then
3013 call readpdb_template(k)
3018 c Distance restraints
3021 C Copy the coordinates from reference coordinates (?)
3025 c write (iout,*) "c(",j,i,") =",c(j,i)
3029 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3031 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3032 open (ientin,file=tpl_k_rescore,status='old')
3033 if (nnt.gt.1) rescore(k,1)=0.0d0
3034 do irec=nnt,nct ! loop for reading res sim
3035 if (read2sigma) then
3036 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3037 & rescore3_tmp,idomain_tmp
3039 idomain(k,i_tmp)=idomain_tmp
3040 rescore(k,i_tmp)=rescore_tmp
3041 rescore2(k,i_tmp)=rescore2_tmp
3042 rescore3(k,i_tmp)=rescore3_tmp
3043 if (.not. out1file .or. me.eq.king)
3044 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3045 & i_tmp,rescore2_tmp,rescore_tmp,
3046 & rescore3_tmp,idomain_tmp
3049 read (ientin,*,end=1401) rescore_tmp
3051 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3052 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3053 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3058 if (waga_dist.ne.0.0d0) then
3066 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3067 c write (iout,*) k,i,j,distal,dist2_cut
3069 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3070 & .and. distal.le.dist2_cut ) then
3076 c write (iout,*) "k",k
3077 c write (iout,*) "i",i," j",j," constr_homology",
3082 if (read2sigma) then
3085 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3087 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3088 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3089 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3091 if (odl(k,ii).le.dist_cut) then
3092 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3095 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3096 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3098 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3099 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3103 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3106 l_homo(k,ii)=.false.
3113 c Theta, dihedral and SC retraints
3115 if (waga_angle.gt.0.0d0) then
3116 c open (ientin,file=tpl_k_sigma_dih,status='old')
3117 c do irec=1,maxres-3 ! loop for reading sigma_dih
3118 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3119 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3120 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3121 c & sigma_dih(k,i+nnt-1)
3126 if (idomain(k,i).eq.0) then
3130 dih(k,i)=phiref(i) ! right?
3131 c read (ientin,*) sigma_dih(k,i) ! original variant
3132 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3133 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3134 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3135 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3136 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3138 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3139 & rescore(k,i-2)+rescore(k,i-3))/4.0
3140 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3141 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3142 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3143 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3144 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3145 c Instead of res sim other local measure of b/b str reliability possible
3146 if (sigma_dih(k,i).ne.0)
3147 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3148 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3153 if (waga_theta.gt.0.0d0) then
3154 c open (ientin,file=tpl_k_sigma_theta,status='old')
3155 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3156 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3157 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3158 c & sigma_theta(k,i+nnt-1)
3163 do i = nnt+2,nct ! right? without parallel.
3164 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3165 c do i=ithet_start,ithet_end ! with FG parallel.
3166 if (idomain(k,i).eq.0) then
3167 sigma_theta(k,i)=0.0
3170 thetatpl(k,i)=thetaref(i)
3171 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3172 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3173 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3174 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3175 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3176 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3177 & rescore(k,i-2))/3.0
3178 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3179 if (sigma_theta(k,i).ne.0)
3180 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3182 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3183 c rescore(k,i-2) ! right expression ?
3184 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3188 if (waga_d.gt.0.0d0) then
3189 c open (ientin,file=tpl_k_sigma_d,status='old')
3190 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3191 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3192 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3193 c & sigma_d(k,i+nnt-1)
3197 do i = nnt,nct ! right? without parallel.
3198 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3199 c do i=loc_start,loc_end ! with FG parallel.
3200 if (itype(i).eq.10) cycle
3201 if (idomain(k,i).eq.0 ) then
3208 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3209 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3210 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3211 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3212 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3213 if (sigma_d(k,i).ne.0)
3214 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3216 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3217 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3218 c read (ientin,*) sigma_d(k,i) ! 1st variant
3223 c remove distance restraints not used in any model from the list
3224 c shift data in all arrays
3226 if (waga_dist.ne.0.0d0) then
3232 if (ii_in_use(ii).eq.0.and.liiflag) then
3236 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3237 & .not.liiflag.and.ii.eq.lim_odl) then
3238 if (ii.eq.lim_odl) then
3239 iishift=ii-iistart+1
3244 do ki=iistart,lim_odl-iishift
3245 ires_homo(ki)=ires_homo(ki+iishift)
3246 jres_homo(ki)=jres_homo(ki+iishift)
3247 ii_in_use(ki)=ii_in_use(ki+iishift)
3248 do k=1,constr_homology
3249 odl(k,ki)=odl(k,ki+iishift)
3250 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3251 l_homo(k,ki)=l_homo(k,ki+iishift)
3255 lim_odl=lim_odl-iishift
3261 endif ! .not. klapaucjusz
3263 if (constr_homology.gt.0) call homology_partition
3264 if (constr_homology.gt.0) call init_int_table
3265 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3266 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3270 if (.not.lprn) return
3271 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3272 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3273 write (iout,*) "Distance restraints from templates"
3275 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3276 & ii,ires_homo(ii),jres_homo(ii),
3277 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3278 & ki=1,constr_homology)
3280 write (iout,*) "Dihedral angle restraints from templates"
3282 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3283 & (rad2deg*dih(ki,i),
3284 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3286 write (iout,*) "Virtual-bond angle restraints from templates"
3288 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3289 & (rad2deg*thetatpl(ki,i),
3290 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3292 write (iout,*) "SC restraints from templates"
3294 write(iout,'(i5,100(4f8.2,4x))') i,
3295 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3296 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3299 c -----------------------------------------------------------------
3302 c----------------------------------------------------------------------
3305 subroutine flush(iu)
3310 subroutine flush(iu)
3315 c------------------------------------------------------------------------------
3316 subroutine copy_to_tmp(source)
3317 include "DIMENSIONS"
3318 include "COMMON.IOUNITS"
3319 character*(*) source
3320 character* 256 tmpfile
3324 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3325 inquire(file=tmpfile,exist=ex)
3327 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3328 & " to temporary directory..."
3329 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3330 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3334 c------------------------------------------------------------------------------
3335 subroutine move_from_tmp(source)
3336 include "DIMENSIONS"
3337 include "COMMON.IOUNITS"
3338 character*(*) source
3341 write (*,*) "Moving ",source(:ilen(source)),
3342 & " from temporary directory to working directory"
3343 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3344 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3347 c------------------------------------------------------------------------------
3348 subroutine random_init(seed)
3350 C Initialize random number generator
3352 implicit real*8 (a-h,o-z)
3353 include 'DIMENSIONS'
3359 logical OKRandom, prng_restart
3361 integer iseed_array(4)
3363 include 'COMMON.IOUNITS'
3364 include 'COMMON.TIME1'
3365 include 'COMMON.THREAD'
3366 include 'COMMON.SBRIDGE'
3367 include 'COMMON.CONTROL'
3368 include 'COMMON.MCM'
3369 include 'COMMON.MAP'
3370 include 'COMMON.HEADER'
3371 include 'COMMON.CSA'
3372 include 'COMMON.CHAIN'
3373 include 'COMMON.MUCA'
3375 include 'COMMON.FFIELD'
3376 include 'COMMON.SETUP'
3377 iseed=-dint(dabs(seed))
3378 if (iseed.eq.0) then
3379 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3380 & 'Random seed undefined. The program will stop.'
3381 write (*,'(/80(1h*)/20x,a/80(1h*))')
3382 & 'Random seed undefined. The program will stop.'
3384 call mpi_finalize(mpi_comm_world,ierr)
3386 stop 'Bad random seed.'
3389 if (fg_rank.eq.0) then
3393 if(me.eq.king .or. .not. out1file)
3394 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3395 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3396 OKRandom = prng_restart(me,iseedi8)
3399 tmp=65536.0d0**(4-i)
3400 iseed_array(i) = dint(seed/tmp)
3401 seed=seed-iseed_array(i)*tmp
3404 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3405 & (iseed_array(i),i=1,4)
3406 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3407 & (iseed_array(i),i=1,4)
3408 OKRandom = prng_restart(me,iseed_array)
3411 c r1 = prng_next(me)
3412 r1=ran_number(0.0D0,1.0D0)
3413 if(me.eq.king .or. .not. out1file)
3414 & write (iout,*) 'ran_num',r1
3415 if (r1.lt.0.0d0) OKRandom=.false.
3417 if (.not.OKRandom) then
3418 write (iout,*) 'PRNG IS NOT WORKING!!!'
3419 print *,'PRNG IS NOT WORKING!!!'
3422 call mpi_abort(mpi_comm_world,error_msg,ierr)
3425 write (iout,*) 'too many processors for parallel prng'
3426 write (*,*) 'too many processors for parallel prng'
3434 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3439 c -----------------------------------------------------------------
3440 subroutine read_klapaucjusz
3442 include 'DIMENSIONS'
3446 include 'COMMON.SETUP'
3447 include 'COMMON.CONTROL'
3448 include 'COMMON.CHAIN'
3449 include 'COMMON.IOUNITS'
3451 include 'COMMON.GEO'
3452 include 'COMMON.INTERACT'
3453 include 'COMMON.NAMES'
3454 character*256 fragfile
3455 integer ninclust(maxclust),inclust(max_template,maxclust),
3456 & nresclust(maxclust),iresclust(maxres,maxclust)
3459 character*24 model_ki_dist, model_ki_angle
3460 character*500 controlcard
3461 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
3462 logical lprn /.true./
3468 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3469 double precision, dimension (max_template,maxres) :: rescore
3470 double precision, dimension (max_template,maxres) :: rescore2
3471 character*24 pdbfile,tpl_k_rescore
3474 c For new homol impl
3476 include 'COMMON.VAR'
3478 call getenv("FRAGFILE",fragfile)
3479 open(ientin,file=fragfile,status="old",err=10)
3480 read(ientin,*) constr_homology,nclust
3486 do k=1,constr_homology
3487 read(ientin,'(a)') pdbfile
3488 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3489 & pdbfile(:ilen(pdbfile))
3490 open(ipdbin,file=pdbfile,status='old',err=33)
3492 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3493 & pdbfile(:ilen(pdbfile))
3497 call readpdb_template(k)
3505 read(ientin,*) ninclust(i),nresclust(i)
3506 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3507 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3510 c Loop over clusters
3513 do ll = 1,ninclust(l)
3521 idomain(k,iresclust(i,l)+1) = 1
3523 idomain(k,iresclust(i,l)) = 1
3527 c Distance restraints
3530 C Copy the coordinates from reference coordinates (?)
3534 c write (iout,*) "c(",j,i,") =",c(j,i)
3537 call int_from_cart(.true.,.false.)
3538 call sc_loc_geom(.false.)
3540 thetaref(i)=theta(i)
3543 if (waga_dist.ne.0.0d0) then
3551 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3552 c write (iout,*) k,i,j,distal,dist2_cut
3554 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3555 & .and. distal.le.dist2_cut ) then
3561 c write (iout,*) "k",k
3562 c write (iout,*) "i",i," j",j," constr_homology",
3567 if (read2sigma) then
3570 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3572 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3573 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3574 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3576 if (odl(k,ii).le.dist_cut) then
3577 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3580 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3581 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3583 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3584 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3588 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3591 c l_homo(k,ii)=.false.
3598 c Theta, dihedral and SC retraints
3600 if (waga_angle.gt.0.0d0) then
3602 if (idomain(k,i).eq.0) then
3603 c sigma_dih(k,i)=0.0
3607 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3608 & rescore(k,i-2)+rescore(k,i-3))/4.0
3609 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3610 c & " sigma_dihed",sigma_dih(k,i)
3611 if (sigma_dih(k,i).ne.0)
3612 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3617 if (waga_theta.gt.0.0d0) then
3619 if (idomain(k,i).eq.0) then
3620 c sigma_theta(k,i)=0.0
3623 thetatpl(k,i)=thetaref(i)
3624 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3625 & rescore(k,i-2))/3.0
3626 if (sigma_theta(k,i).ne.0)
3627 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3631 if (waga_d.gt.0.0d0) then
3633 if (itype(i).eq.10) cycle
3634 if (idomain(k,i).eq.0 ) then
3641 sigma_d(k,i)=rescore(k,i)
3642 if (sigma_d(k,i).ne.0)
3643 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3644 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3650 c remove distance restraints not used in any model from the list
3651 c shift data in all arrays
3653 if (waga_dist.ne.0.0d0) then
3659 if (ii_in_use(ii).eq.0.and.liiflag) then
3663 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3664 & .not.liiflag.and.ii.eq.lim_odl) then
3665 if (ii.eq.lim_odl) then
3666 iishift=ii-iistart+1
3671 do ki=iistart,lim_odl-iishift
3672 ires_homo(ki)=ires_homo(ki+iishift)
3673 jres_homo(ki)=jres_homo(ki+iishift)
3674 ii_in_use(ki)=ii_in_use(ki+iishift)
3675 do k=1,constr_homology
3676 odl(k,ki)=odl(k,ki+iishift)
3677 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3678 l_homo(k,ki)=l_homo(k,ki+iishift)
3682 lim_odl=lim_odl-iishift
3689 10 stop "Error infragment file"
3691 c----------------------------------------------------------------------