8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
14 C Read job setup parameters
16 C Read force-field parameters except weights
18 C Read control parameters for energy minimzation if required
19 if (minim) call read_minim
20 C Read MCM control parameters if required
21 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
22 C Read MD control parameters if reqjuired
23 if (modecalc.eq.12) call read_MDpar
24 C Read MREMD control parameters if required
25 if (modecalc.eq.14) then
29 C Read MUCA control parameters if required
30 if (lmuca) call read_muca
31 C Read CSA control parameters if required (from fort.40 if exists
32 C otherwise from general input file)
33 if (modecalc.eq.8) then
34 inquire (file="fort.40",exist=file_exist)
35 if (.not.file_exist) call csaread
37 cfmc if (modecalc.eq.10) call mcmfread
38 c Read energy-term weights and disulfide parameters
40 C Read molecule information, molecule geometry, energy-term weights, and
41 C restraints if requested
43 C Print restraint information
45 if (.not. out1file .or. me.eq.king) then
48 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
49 & "The following",nhpb-nss,
50 & " distance restraints have been imposed:",
51 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
54 write (iout,'(4i5,2f8.2,3f10.5,2i5)')i-nss,ihpb(i),jhpb(i),
55 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i),fordepth(i),xlscore(i),
60 write (iout,'(/a,i5,a/4a5,2a8,3a10,2a5)')
61 & "The following",npeak,
62 & " NMR peak restraints have been imposed:",
63 & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V",
64 & " score"," type"," ipeak"
66 do j=ipeak(1,i),ipeak(2,i)
67 write (iout,'(5i5,2f8.2,2f10.5,i5)')i,j,ihpb_peak(j),
68 & jhpb_peak(j),ibecarb_peak(j),dhpb_peak(j),dhpb1_peak(j),
69 & forcon_peak(j),fordepth_peak(i),irestr_type_peak(j)
72 write (iout,*) "The ipeak array"
74 write (iout,'(3i5)' ) i,ipeak(1,i),ipeak(2,i)
80 c print *,"Processor",myrank," leaves READRTNS"
83 C-------------------------------------------------------------------------------
84 subroutine read_control
92 logical OKRandom, prng_restart
95 include 'COMMON.IOUNITS'
96 include 'COMMON.TIME1'
97 include 'COMMON.THREAD'
98 include 'COMMON.SBRIDGE'
99 include 'COMMON.CONTROL'
100 include 'COMMON.SAXS'
103 include 'COMMON.HEADER'
105 include 'COMMON.CHAIN'
106 include 'COMMON.MUCA'
108 include 'COMMON.FFIELD'
109 include 'COMMON.INTERACT'
110 include 'COMMON.SETUP'
111 include 'COMMON.SPLITELE'
112 include 'COMMON.SHIELD'
115 integer KDIAG,ICORFL,IXDR
116 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
117 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
119 character*320 controlcard
120 double precision seed
125 read (INP,'(a)') titel
126 call card_concat(controlcard)
127 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
128 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
129 call reada(controlcard,'SEED',seed,0.0D0)
130 call random_init(seed)
131 C Set up the time limit (caution! The time must be input in minutes!)
132 read_cart=index(controlcard,'READ_CART').gt.0
133 out_cart=index(controlcard,'OUT_CART').gt.0
134 out_int=index(controlcard,'OUT_INT').gt.0
135 gmatout=index(controlcard,'GMATOUT').gt.0
136 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
137 C this variable with_theta_constr is the variable which allow to read and execute the
138 C constrains on theta angles WITH_THETA_CONSTR is the keyword
139 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
140 write (iout,*) "with_theta_constr ",with_theta_constr
141 call readi(controlcard,'NSAXS',nsaxs,0)
142 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
143 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
144 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
145 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
146 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
147 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
148 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
149 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
150 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
151 call readi(controlcard,'SYM',symetr,1)
152 call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours
153 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
154 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
155 c call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
156 c call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
157 c call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
158 c call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
159 c call reada(controlcard,'DRMS',drms,0.1D0)
160 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
161 c write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
162 c write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
163 c write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
164 c write (iout,'(a,f10.1)')'DRMS = ',drms
165 cc write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
166 c write (iout,'(a,f10.1)') 'Time limit (min):',timlim
168 call readi(controlcard,'NZ_START',nz_start,0)
169 call readi(controlcard,'NZ_END',nz_end,0)
170 c call readi(controlcard,'IZ_SC',iz_sc,0)
172 safety = 60.0d0*safety
174 call reada(controlcard,"T_BATH",t_bath,300.0d0)
175 minim=(index(controlcard,'MINIMIZE').gt.0)
176 dccart=(index(controlcard,'CART').gt.0)
177 overlapsc=(index(controlcard,'OVERLAP').gt.0)
178 overlapsc=.not.overlapsc
179 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
180 searchsc=.not.searchsc
181 sideadd=(index(controlcard,'SIDEADD').gt.0)
182 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
183 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
184 outpdb=(index(controlcard,'PDBOUT').gt.0)
185 outmol2=(index(controlcard,'MOL2OUT').gt.0)
186 pdbref=(index(controlcard,'PDBREF').gt.0)
187 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
188 indpdb=index(controlcard,'PDBSTART')
189 extconf=(index(controlcard,'EXTCONF').gt.0)
190 AFMlog=(index(controlcard,'AFM'))
191 selfguide=(index(controlcard,'SELFGUIDE'))
192 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
193 call readi(controlcard,'TUBEMOD',tubelog,0)
194 c write (iout,*) TUBElog,"TUBEMODE"
195 call readi(controlcard,'IPRINT',iprint,0)
196 C SHIELD keyword sets if the shielding effect of side-chains is used
197 C 0 denots no shielding is used all peptide are equally despite the
198 C solvent accesible area
199 C 1 the newly introduced function
200 C 2 reseved for further possible developement
201 call readi(controlcard,'SHIELD',shield_mode,0)
202 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
203 write(iout,*) "shield_mode",shield_mode
205 call readi(controlcard,'TORMODE',tor_mode,0)
206 C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
207 write(iout,*) "torsional and valence angle mode",tor_mode
208 call readi(controlcard,'MAXGEN',maxgen,10000)
209 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
210 call readi(controlcard,"KDIAG",kdiag,0)
211 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
212 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
213 & write (iout,*) "RESCALE_MODE",rescale_mode
214 split_ene=index(controlcard,'SPLIT_ENE').gt.0
215 if (index(controlcard,'REGULAR').gt.0.0D0) then
216 call reada(controlcard,'WEIDIS',weidis,0.1D0)
221 if (index(controlcard,"CASC").gt.0) then
223 c else if (index(controlcard,"CAONLY").gt.0) then
225 else if (index(controlcard,"SCONLY").gt.0) then
229 c write (iout,*) "ERROR! Must specify CASC CAONLY or SCONLY when",
230 c & " specifying REFSTR, PDBREF or REGULAR."
234 if (index(controlcard,'CHECKGRAD').gt.0) then
236 if (index(controlcard,' CART').gt.0) then
238 elseif (index(controlcard,'CARINT').gt.0) then
243 call reada(controlcard,'DELTA',aincr,1.0d-7)
244 c write (iout,*) "icheckgrad",icheckgrad
245 elseif (index(controlcard,'THREAD').gt.0) then
247 call readi(controlcard,'THREAD',nthread,0)
248 if (nthread.gt.0) then
249 call reada(controlcard,'WEIDIS',weidis,0.1D0)
252 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
253 stop 'Error termination in Read_Control.'
255 else if (index(controlcard,'MCMA').gt.0) then
257 else if (index(controlcard,'MCEE').gt.0) then
259 else if (index(controlcard,'MULTCONF').gt.0) then
261 else if (index(controlcard,'MAP').gt.0) then
263 call readi(controlcard,'MAP',nmap,0)
264 else if (index(controlcard,'CSA').gt.0) then
266 crc else if (index(controlcard,'ZSCORE').gt.0) then
268 crc ZSCORE is rm from UNRES, modecalc=9 is available
271 cfcm else if (index(controlcard,'MCMF').gt.0) then
273 else if (index(controlcard,'SOFTREG').gt.0) then
275 else if (index(controlcard,'CHECK_BOND').gt.0) then
277 else if (index(controlcard,'TEST').gt.0) then
279 else if (index(controlcard,'MD').gt.0) then
281 else if (index(controlcard,'RE ').gt.0) then
285 lmuca=index(controlcard,'MUCA').gt.0
286 call readi(controlcard,'MUCADYN',mucadyn,0)
287 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
288 if (lmuca .and. (me.eq.king .or. .not.out1file ))
290 write (iout,*) 'MUCADYN=',mucadyn
291 write (iout,*) 'MUCASMOOTH=',muca_smooth
294 iscode=index(controlcard,'ONE_LETTER')
295 indphi=index(controlcard,'PHI')
296 indback=index(controlcard,'BACK')
297 iranconf=index(controlcard,'RAND_CONF')
298 i2ndstr=index(controlcard,'USE_SEC_PRED')
299 gradout=index(controlcard,'GRADOUT').gt.0
300 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
301 C DISTCHAINMAX become obsolete for periodic boundry condition
302 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
303 C Reading the dimensions of box in x,y,z coordinates
304 call reada(controlcard,'BOXX',boxxsize,100.0d0)
305 call reada(controlcard,'BOXY',boxysize,100.0d0)
306 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
307 c Cutoff range for interactions
308 call reada(controlcard,"R_CUT",r_cut,15.0d0)
309 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
310 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
311 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
312 if (lipthick.gt.0.0d0) then
313 bordliptop=(boxzsize+lipthick)/2.0
314 bordlipbot=bordliptop-lipthick
316 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
317 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
318 buflipbot=bordlipbot+lipbufthick
319 bufliptop=bordliptop-lipbufthick
320 if ((lipbufthick*2.0d0).gt.lipthick)
321 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
323 write(iout,*) "bordliptop=",bordliptop
324 write(iout,*) "bordlipbot=",bordlipbot
325 write(iout,*) "bufliptop=",bufliptop
326 write(iout,*) "buflipbot=",buflipbot
327 write (iout,*) "SHIELD MODE",shield_mode
328 if (TUBElog.gt.0) then
329 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
330 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
331 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
332 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
333 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
334 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
335 buftubebot=bordtubebot+tubebufthick
336 buftubetop=bordtubetop-tubebufthick
338 if (me.eq.king .or. .not.out1file )
339 & write (iout,*) "DISTCHAINMAX",distchainmax
341 if(me.eq.king.or..not.out1file)
342 & write (iout,'(2a)') diagmeth(kdiag),
343 & ' routine used to diagonalize matrices.'
346 c--------------------------------------------------------------------------
347 subroutine read_REMDpar
353 include 'COMMON.IOUNITS'
354 include 'COMMON.TIME1'
357 include 'COMMON.LANGEVIN'
360 include 'COMMON.LANGEVIN.lang0.5diag'
362 include 'COMMON.LANGEVIN.lang0'
365 include 'COMMON.INTERACT'
366 include 'COMMON.NAMES'
368 include 'COMMON.REMD'
369 include 'COMMON.CONTROL'
370 include 'COMMON.SETUP'
372 character*320 controlcard
373 character*3200 controlcard1
374 integer iremd_m_total,i
376 if(me.eq.king.or..not.out1file)
377 & write (iout,*) "REMD setup"
379 call card_concat(controlcard)
380 call readi(controlcard,"NREP",nrep,3)
381 call readi(controlcard,"NSTEX",nstex,1000)
382 call reada(controlcard,"RETMIN",retmin,10.0d0)
383 call reada(controlcard,"RETMAX",retmax,1000.0d0)
384 mremdsync=(index(controlcard,'SYNC').gt.0)
385 call readi(controlcard,"NSYN",i_sync_step,100)
386 restart1file=(index(controlcard,'REST1FILE').gt.0)
387 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
388 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
389 if(max_cache_traj_use.gt.max_cache_traj)
390 & max_cache_traj_use=max_cache_traj
391 if(me.eq.king.or..not.out1file) then
392 cd if (traj1file) then
393 crc caching is in testing - NTWX is not ignored
394 cd write (iout,*) "NTWX value is ignored"
395 cd write (iout,*) " trajectory is stored to one file by master"
396 cd write (iout,*) " before exchange at NSTEX intervals"
398 write (iout,*) "NREP= ",nrep
399 write (iout,*) "NSTEX= ",nstex
400 write (iout,*) "SYNC= ",mremdsync
401 write (iout,*) "NSYN= ",i_sync_step
402 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
405 if (index(controlcard,'TLIST').gt.0) then
407 call card_concat(controlcard1)
408 read(controlcard1,*) (remd_t(i),i=1,nrep)
409 if(me.eq.king.or..not.out1file)
410 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
413 if (index(controlcard,'MLIST').gt.0) then
415 call card_concat(controlcard1)
416 read(controlcard1,*) (remd_m(i),i=1,nrep)
417 if(me.eq.king.or..not.out1file) then
418 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
421 iremd_m_total=iremd_m_total+remd_m(i)
423 write (iout,*) 'Total number of replicas ',iremd_m_total
428 & "Adaptive (PMF-biased) umbrella sampling will be run"
431 if(me.eq.king.or..not.out1file)
432 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
435 c--------------------------------------------------------------------------
436 subroutine read_MDpar
442 include 'COMMON.IOUNITS'
443 include 'COMMON.TIME1'
445 include 'COMMON.QRESTR'
447 include 'COMMON.LANGEVIN'
450 include 'COMMON.LANGEVIN.lang0.5diag'
452 include 'COMMON.LANGEVIN.lang0'
455 include 'COMMON.INTERACT'
456 include 'COMMON.NAMES'
458 include 'COMMON.SETUP'
459 include 'COMMON.CONTROL'
460 include 'COMMON.SPLITELE'
461 include 'COMMON.FFIELD'
463 character*320 controlcard
467 call card_concat(controlcard)
468 call readi(controlcard,"NSTEP",n_timestep,1000000)
469 call readi(controlcard,"NTWE",ntwe,100)
470 call readi(controlcard,"NTWX",ntwx,1000)
471 call reada(controlcard,"DT",d_time,1.0d-1)
472 call reada(controlcard,"DVMAX",dvmax,2.0d1)
473 call reada(controlcard,"DAMAX",damax,1.0d1)
474 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
475 call readi(controlcard,"LANG",lang,0)
476 RESPA = index(controlcard,"RESPA") .gt. 0
477 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
478 ntime_split0=ntime_split
479 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
480 ntime_split0=ntime_split
481 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
482 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
483 rest = index(controlcard,"REST").gt.0
484 tbf = index(controlcard,"TBF").gt.0
485 usampl = index(controlcard,"USAMPL").gt.0
486 scale_umb = index(controlcard,"SCALE_UMB").gt.0
487 adaptive = index(controlcard,"ADAPTIVE").gt.0
488 mdpdb = index(controlcard,"MDPDB").gt.0
489 call reada(controlcard,"T_BATH",t_bath,300.0d0)
490 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
491 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
492 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
493 if (count_reset_moment.eq.0) count_reset_moment=1000000000
494 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
495 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
496 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
497 if (count_reset_vel.eq.0) count_reset_vel=1000000000
498 large = index(controlcard,"LARGE").gt.0
499 print_compon = index(controlcard,"PRINT_COMPON").gt.0
500 rattle = index(controlcard,"RATTLE").gt.0
501 preminim = index(controlcard,"PREMINIM").gt.0
502 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
506 write (iout,*) "PREMINIM ",preminim
508 if (index(controlcard,'CART').gt.0) dccart=.true.
509 write (iout,*) "dccart ",dccart
510 write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
511 if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
513 c if performing umbrella sampling, fragments constrained are read from the fragment file
516 write (iout,*) "Umbrella sampling will be run"
517 if (scale_umb.and.adaptive) then
518 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
519 write (iout,*) "Select one of those and re-run the job."
522 if (scale_umb) write (iout,*)
523 &"Umbrella-restraint force constants will be scaled by temperature"
527 if(me.eq.king.or..not.out1file) then
529 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
531 write (iout,'(a)') "The units are:"
532 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
533 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
534 & " acceleration: angstrom/(48.9 fs)**2"
535 write (iout,'(a)') "energy: kcal/mol, temperature: K"
537 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
538 write (iout,'(a60,f10.5,a)')
539 & "Initial time step of numerical integration:",d_time,
541 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
543 write (iout,'(2a,i4,a)')
544 & "A-MTS algorithm used; initial time step for fast-varying",
545 & " short-range forces split into",ntime_split," steps."
546 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
547 & r_cut," lambda",rlamb
549 write (iout,'(2a,f10.5)')
550 & "Maximum acceleration threshold to reduce the time step",
551 & "/increase split number:",damax
552 write (iout,'(2a,f10.5)')
553 & "Maximum predicted energy drift to reduce the timestep",
554 & "/increase split number:",edriftmax
555 write (iout,'(a60,f10.5)')
556 & "Maximum velocity threshold to reduce velocities:",dvmax
557 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
558 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
559 if (rattle) write (iout,'(a60)')
560 & "Rattle algorithm used to constrain the virtual bonds"
564 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
565 call reada(controlcard,"RWAT",rwat,1.4d0)
566 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
567 surfarea=index(controlcard,"SURFAREA").gt.0
568 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
569 if(me.eq.king.or..not.out1file)then
570 write (iout,'(/a,$)') "Langevin dynamics calculation"
573 & " with direct integration of Langevin equations"
574 else if (lang.eq.2) then
575 write (iout,'(a/)') " with TINKER stochasic MD integrator"
576 else if (lang.eq.3) then
577 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
578 else if (lang.eq.4) then
579 write (iout,'(a/)') " in overdamped mode"
581 write (iout,'(//a,i5)')
582 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
585 write (iout,'(a60,f10.5)') "Temperature:",t_bath
586 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
587 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
588 write (iout,'(a60,f10.5)')
589 & "Scaling factor of the friction forces:",scal_fric
590 if (surfarea) write (iout,'(2a,i10,a)')
591 & "Friction coefficients will be scaled by solvent-accessible",
592 & " surface area every",reset_fricmat," steps."
594 c Calculate friction coefficients and bounds of stochastic forces
595 eta=6*pi*cPoise*etawat
596 if(me.eq.king.or..not.out1file)
597 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
599 gamp=scal_fric*(pstok+rwat)*eta
600 stdfp=dsqrt(2*Rb*t_bath/d_time)
602 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
603 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
605 if(me.eq.king.or..not.out1file)then
606 write (iout,'(/2a/)')
607 & "Radii of site types and friction coefficients and std's of",
608 & " stochastic forces of fully exposed sites"
609 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
611 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
612 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
616 if(me.eq.king.or..not.out1file)then
617 write (iout,'(a)') "Berendsen bath calculation"
618 write (iout,'(a60,f10.5)') "Temperature:",t_bath
619 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
621 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
622 & count_reset_moment," steps"
624 & write (iout,'(a,i10,a)')
625 & "Velocities will be reset at random every",count_reset_vel,
629 if(me.eq.king.or..not.out1file)
630 & write (iout,'(a31)') "Microcanonical mode calculation"
632 if(me.eq.king.or..not.out1file)then
633 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
635 write(iout,*) "MD running with constraints."
636 write(iout,*) "Equilibration time ", eq_time, " mtus."
637 write(iout,*) "Constraining ", nfrag," fragments."
638 write(iout,*) "Length of each fragment, weight and q0:"
640 write (iout,*) "Set of restraints #",iset
642 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
643 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
645 write(iout,*) "constraints between ", npair, "fragments."
646 write(iout,*) "constraint pairs, weights and q0:"
648 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
649 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
651 write(iout,*) "angle constraints within ", nfrag_back,
652 & "backbone fragments."
654 write(iout,*) "fragment, weights, q0:"
656 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
657 & ifrag_back(2,i,iset),
658 & wfrag_back(1,i,iset),qin_back(1,i,iset),
659 & wfrag_back(2,i,iset),qin_back(2,i,iset),
660 & wfrag_back(3,i,iset),qin_back(3,i,iset)
663 write(iout,*) "fragment, weights:"
665 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
666 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
667 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
671 iset=mod(kolor,nset)+1
674 if(me.eq.king.or..not.out1file)
675 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
678 c------------------------------------------------------------------------------
681 C Read molecular data.
687 integer error_msg,ierror,ierr,ierrcode
689 include 'COMMON.IOUNITS'
692 include 'COMMON.INTERACT'
693 include 'COMMON.LOCAL'
694 include 'COMMON.NAMES'
695 include 'COMMON.CHAIN'
696 include 'COMMON.FFIELD'
697 include 'COMMON.SBRIDGE'
698 include 'COMMON.HEADER'
699 include 'COMMON.CONTROL'
700 include 'COMMON.SAXS'
701 include 'COMMON.DBASE'
702 include 'COMMON.THREAD'
703 include 'COMMON.CONTACTS'
704 include 'COMMON.TORCNSTR'
705 include 'COMMON.TIME1'
706 include 'COMMON.BOUNDS'
708 include 'COMMON.SETUP'
709 include 'COMMON.SHIELD'
710 character*4 sequence(maxres)
712 double precision x(maxvar)
713 character*256 pdbfile
714 character*400 weightcard
715 character*80 weightcard_t,ucase
716 integer itype_pdb(maxres)
717 common /pizda/ itype_pdb
718 logical seq_comp,fail
719 double precision energia(0:n_ene)
720 double precision secprob(3,maxdih_constr)
722 double precision phihel,phibet,sigmahel,sigmabet
723 integer iti,nsi,maxsi
727 integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2
728 double precision sumv
730 C Read PDB structure if applicable
732 if (indpdb.gt.0 .or. pdbref) then
733 read(inp,'(a)') pdbfile
734 if(me.eq.king.or..not.out1file)
735 & write (iout,'(2a)') 'PDB data will be read from file ',
736 & pdbfile(:ilen(pdbfile))
737 open(ipdbin,file=pdbfile,status='old',err=33)
739 33 write (iout,'(a)') 'Error opening PDB file.'
750 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
751 & (crefjlee(j,i+nres),j=1,3)
754 if(me.eq.king.or..not.out1file)
755 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
756 & ' nstart_sup=',nstart_sup
758 itype_pdb(i)=itype(i)
762 nct=nstart_sup+nsup-1
763 call contact(.false.,ncont_ref,icont_ref,co)
766 if(me.eq.king.or..not.out1file)
767 & write(iout,*)'Adding sidechains'
771 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
774 do while (fail.and.nsi.le.maxsi)
775 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
778 if(fail) write(iout,*)'Adding sidechain failed for res ',
779 & i,' after ',nsi,' trials'
784 if (indpdb.eq.0) then
785 C Read sequence if not taken from the pdb file.
787 c print *,'nres=',nres
788 if (iscode.gt.0) then
789 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
791 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
793 C Convert sequence to numeric code
795 itype(i)=rescode(i,sequence(i),iscode)
799 c print '(20i4)',(itype(i),i=1,nres)
802 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
804 if (itype(i).eq.ntyp1) then
808 else if (iabs(itype(i+1)).ne.20) then
810 else if (iabs(itype(i)).ne.20) then
817 if(me.eq.king.or..not.out1file)then
818 write (iout,*) "ITEL"
820 write (iout,*) i,itype(i),itel(i)
822 print *,'Call Read_Bridge.'
826 cd print *,'NNT=',NNT,' NCT=',NCT
827 call seq2chains(nres,itype,nchain,chain_length,chain_border,
830 chain_border1(2,1)=chain_border(2,1)+1
832 chain_border1(1,i)=chain_border(1,i)-1
833 chain_border1(2,i)=chain_border(2,i)+1
835 chain_border1(1,nchain)=chain_border(1,nchain)-1
836 chain_border1(2,nchain)=nres
837 write(iout,*) "nres",nres," nchain",nchain
839 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
840 & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
842 call chain_symmetry(nchain,nres,itype,chain_border,
843 & chain_length,npermchain,tabpermchain)
845 c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
848 write(iout,*) "residue permutations"
850 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
853 if (itype(1).eq.ntyp1) nnt=2
854 if (itype(nres).eq.ntyp1) nct=nct-1
855 write (iout,*) "nnt",nnt," nct",nct
858 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
859 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
861 print*, 'init_dfa_vars finished!'
863 print*, 'read_dfa_info finished!'
867 if(me.eq.king.or..not.out1file)
868 & write (iout,'(a,i3)') 'nsup=',nsup
870 if (nsup.le.(nct-nnt+1)) then
871 do i=0,nct-nnt+1-nsup
872 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
878 & 'Error - sequences to be superposed do not match.'
881 do i=0,nsup-(nct-nnt+1)
882 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
884 nstart_sup=nstart_sup+i
890 & 'Error - sequences to be superposed do not match.'
893 if (nsup.eq.0) nsup=nct-nnt
894 if (nstart_sup.eq.0) nstart_sup=nnt
895 if (nstart_seq.eq.0) nstart_seq=nnt
896 if(me.eq.king.or..not.out1file)
897 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
898 & ' nstart_seq=',nstart_seq
901 C 8/13/98 Set limits to generating the dihedral angles
906 read (inp,*) ndih_constr
907 write (iout,*) "ndih_constr",ndih_constr
908 if (ndih_constr.gt.0) then
911 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
914 if(me.eq.king.or..not.out1file)then
917 & 'There are',ndih_constr,' restraints on gamma angles.'
919 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
926 phi0(i)=deg2rad*phi0(i)
927 drange(i)=deg2rad*drange(i)
929 C if(me.eq.king.or..not.out1file)
930 C & write (iout,*) 'FTORS',ftors
933 phibound(1,ii) = phi0(i)-drange(i)
934 phibound(2,ii) = phi0(i)+drange(i)
937 if (me.eq.king .or. .not.out1file) then
939 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
941 write (iout,'(a3,i5,2f10.1)')
942 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
947 else if (ndih_constr.lt.0) then
949 call card_concat(weightcard)
950 call reada(weightcard,"PHIHEL",phihel,50.0D0)
951 call reada(weightcard,"PHIBET",phibet,180.0D0)
952 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
953 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
954 call reada(weightcard,"WDIHC",wdihc,0.591D0)
955 write (iout,*) "Weight of dihedral angle restraints",wdihc
956 read(inp,'(9x,3f7.3)')
957 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
958 write (iout,*) "The secprob array"
960 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
964 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
965 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
966 ndih_constr=ndih_constr+1
967 idih_constr(ndih_constr)=i
970 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
971 sumv=sumv+vpsipred(j,ndih_constr)
974 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
976 phibound(1,ndih_constr)=phihel*deg2rad
977 phibound(2,ndih_constr)=phibet*deg2rad
978 sdihed(1,ndih_constr)=sigmahel*deg2rad
979 sdihed(2,ndih_constr)=sigmabet*deg2rad
983 if(me.eq.king.or..not.out1file)then
986 & 'There are',ndih_constr,
987 & ' bimodal restraints on gamma angles.'
989 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
990 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
991 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
992 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
993 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
994 & (vpsipred(j,i),j=1,3)
1001 C first setting the theta boundaries to 0 to pi
1002 C this mean that there is no energy penalty for any angle occuring this can be applied
1003 C for generate random conformation but is not implemented in this way
1006 C thetabound(2,i)=pi
1008 C begin reading theta constrains this is quartic constrains allowing to
1009 C have smooth second derivative
1010 if (with_theta_constr) then
1011 C with_theta_constr is keyword allowing for occurance of theta constrains
1012 read (inp,*) ntheta_constr
1013 C ntheta_constr is the number of theta constrains
1014 if (ntheta_constr.gt.0) then
1015 C read (inp,*) ftors
1016 read (inp,*) (itheta_constr(i),theta_constr0(i),
1017 & theta_drange(i),for_thet_constr(i),
1018 & i=1,ntheta_constr)
1019 C the above code reads from 1 to ntheta_constr
1020 C itheta_constr(i) residue i for which is theta_constr
1021 C theta_constr0 the global minimum value
1022 C theta_drange is range for which there is no energy penalty
1023 C for_thet_constr is the force constant for quartic energy penalty
1025 if(me.eq.king.or..not.out1file)then
1027 & 'There are',ntheta_constr,' constraints on phi angles.'
1028 do i=1,ntheta_constr
1029 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1031 & for_thet_constr(i)
1034 do i=1,ntheta_constr
1035 theta_constr0(i)=deg2rad*theta_constr0(i)
1036 theta_drange(i)=deg2rad*theta_drange(i)
1038 C if(me.eq.king.or..not.out1file)
1039 C & write (iout,*) 'FTORS',ftors
1040 C do i=1,ntheta_constr
1041 C ii = itheta_constr(i)
1042 C thetabound(1,ii) = phi0(i)-drange(i)
1043 C thetabound(2,ii) = phi0(i)+drange(i)
1045 endif ! ntheta_constr.gt.0
1046 endif! with_theta_constr
1048 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1049 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1050 c--- Zscore rms -------
1051 if (nz_start.eq.0) nz_start=nnt
1052 if (nz_end.eq.0 .and. nsup.gt.0) then
1054 else if (nz_end.eq.0) then
1057 if(me.eq.king.or..not.out1file)then
1058 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1059 write (iout,*) 'IZ_SC=',iz_sc
1061 c----------------------
1064 if (.not.pdbref) then
1065 call read_angles(inp,*38)
1068 38 write (iout,'(a)') 'Error reading reference structure.'
1070 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1071 stop 'Error reading reference structure'
1073 39 call chainbuild_extconf
1075 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1084 call contact(.true.,ncont_ref,icont_ref,co)
1088 write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1090 if (constr_dist.gt.0) call read_dist_constr
1091 write (iout,*) "After read_dist_constr nhpb",nhpb
1092 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1094 call NMRpeak_partition
1095 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1097 if(me.eq.king.or..not.out1file)
1098 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1101 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1103 if(me.eq.king.or..not.out1file)
1104 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1105 & icont_ref(1,i),' ',
1106 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1109 write (iout,*) "calling read_saxs_consrtr",nsaxs
1110 if (nsaxs.gt.0) call read_saxs_constr
1112 if (constr_homology.gt.0) then
1113 call read_constr_homology
1114 if (indpdb.gt.0 .or. pdbref) then
1117 c(j,i)=crefjlee(j,i)
1118 cref(j,i)=crefjlee(j,i)
1123 write (iout,*) "sc_loc_geom: Array C"
1125 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1126 & (c(j,i+nres),j=1,3)
1128 write (iout,*) "Array Cref"
1130 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1131 & (cref(j,i+nres),j=1,3)
1134 call int_from_cart1(.false.)
1135 call sc_loc_geom(.false.)
1137 thetaref(i)=theta(i)
1142 dc(j,i)=c(j,i+1)-c(j,i)
1143 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1148 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1149 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1158 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1159 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1160 & modecalc.ne.10) then
1161 C If input structure hasn't been supplied from the PDB file read or generate
1163 if (iranconf.eq.0 .and. .not. extconf) then
1164 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1165 & write (iout,'(a)') 'Initial geometry will be read in.'
1167 read(inp,'(8f10.5)',end=36,err=36)
1168 & ((c(l,k),l=1,3),k=1,nres),
1169 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1170 write (iout,*) "Exit READ_CART"
1171 c write (iout,'(8f10.5)')
1172 c & ((c(l,k),l=1,3),k=1,nres),
1173 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1177 dc(j,i)=c(j,i+1)-c(j,i)
1178 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1182 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1184 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1185 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1190 c dc_norm(j,i+nres)=0.0d0
1194 call int_from_cart1(.true.)
1195 write (iout,*) "Finish INT_TO_CART"
1196 c write (iout,*) "DC"
1198 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1199 c & (dc(j,i+nres),j=1,3)
1205 call read_angles(inp,*36)
1207 call chainbuild_extconf
1210 36 write (iout,'(a)') 'Error reading angle file.'
1212 call mpi_finalize( MPI_COMM_WORLD,IERR )
1214 stop 'Error reading angle file.'
1216 else if (extconf) then
1217 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1218 & write (iout,'(a)') 'Extended chain initial geometry.'
1220 theta(i)=90d0*deg2rad
1223 phi(i)=180d0*deg2rad
1226 alph(i)=110d0*deg2rad
1229 omeg(i)=-120d0*deg2rad
1230 if (itype(i).le.0) omeg(i)=-omeg(i)
1233 call chainbuild_extconf
1235 if(me.eq.king.or..not.out1file)
1236 & write (iout,'(a)') 'Random-generated initial geometry.'
1239 if (me.eq.king .or. fg_rank.eq.0 .and. (
1240 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1244 call gen_rand_conf(itmp,*30)
1246 30 write (iout,*) 'Failed to generate random conformation',
1247 & ', itrial=',itrial
1248 write (*,*) 'Processor:',me,
1249 & ' Failed to generate random conformation',
1259 write (iout,'(a,i3,a)') 'Processor:',me,
1260 & ' error in generating random conformation.'
1261 write (*,'(a,i3,a)') 'Processor:',me,
1262 & ' error in generating random conformation.'
1265 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1270 & ' error in generating random conformation.'
1275 elseif (modecalc.eq.4) then
1276 read (inp,'(a)') intinname
1277 open (intin,file=intinname,status='old',err=333)
1278 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1279 & write (iout,'(a)') 'intinname',intinname
1280 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1282 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1284 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1286 stop 'Error opening angle file.'
1290 C Generate distance constraints, if the PDB structure is to be regularized.
1291 if (nthread.gt.0) then
1292 call read_threadbase
1295 if (me.eq.king .or. .not. out1file)
1297 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1298 write (iout,'(/a,i3,a)')
1299 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1300 write (iout,'(20i4)') (iss(i),i=1,ns)
1302 write(iout,*)"Running with dynamic disulfide-bond formation"
1304 write (iout,'(/a/)') 'Pre-formed links are:'
1310 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1311 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1317 if (ns.gt.0.and.dyn_ss) then
1321 forcon(i-nss)=forcon(i)
1328 dyn_ss_mask(iss(i))=.true.
1333 c write (iout,*) "DC"
1335 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1336 c & (dc(j,i+nres),j=1,3)
1338 if (i2ndstr.gt.0) call secstrp2dihc
1339 c call geom_to_var(nvar,x)
1340 c call etotal(energia(0))
1341 c call enerprint(energia(0))
1342 c call briefout(0,etot)
1344 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1345 cd write (iout,'(a)') 'Variable list:'
1346 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1348 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1349 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1350 & 'Processor',myrank,': end reading molecular data.'
1355 c--------------------------------------------------------------------------
1356 logical function seq_comp(itypea,itypeb,length)
1358 integer length,itypea(length),itypeb(length)
1361 if (itypea(i).ne.itypeb(i)) then
1369 c-----------------------------------------------------------------------------
1370 subroutine read_bridge
1371 C Read information about disulfide bridges.
1373 include 'DIMENSIONS'
1378 include 'COMMON.IOUNITS'
1379 include 'COMMON.GEO'
1380 include 'COMMON.VAR'
1381 include 'COMMON.INTERACT'
1382 include 'COMMON.LOCAL'
1383 include 'COMMON.NAMES'
1384 include 'COMMON.CHAIN'
1385 include 'COMMON.FFIELD'
1386 include 'COMMON.SBRIDGE'
1387 include 'COMMON.HEADER'
1388 include 'COMMON.CONTROL'
1389 include 'COMMON.DBASE'
1390 include 'COMMON.THREAD'
1391 include 'COMMON.TIME1'
1392 include 'COMMON.SETUP'
1394 C Read bridging residues.
1395 read (inp,*) ns,(iss(i),i=1,ns)
1397 if(me.eq.king.or..not.out1file)
1398 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1399 C Check whether the specified bridging residues are cystines.
1401 if (itype(iss(i)).ne.1) then
1402 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1403 & 'Do you REALLY think that the residue ',
1404 & restyp(itype(iss(i))),i,
1405 & ' can form a disulfide bridge?!!!'
1406 write (*,'(2a,i3,a)')
1407 & 'Do you REALLY think that the residue ',
1408 & restyp(itype(iss(i))),i,
1409 & ' can form a disulfide bridge?!!!'
1411 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1416 C Read preformed bridges.
1418 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1420 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1423 C Check if the residues involved in bridges are in the specified list of
1424 C bridging residues.
1427 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1428 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1429 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1430 & ' contains residues present in other pairs.'
1431 write (*,'(a,i3,a)') 'Disulfide pair',i,
1432 & ' contains residues present in other pairs.'
1434 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1440 if (ihpb(i).eq.iss(j)) goto 10
1442 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1445 if (jhpb(i).eq.iss(j)) goto 20
1447 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1453 ihpb(i)=ihpb(i)+nres
1454 jhpb(i)=jhpb(i)+nres
1460 c----------------------------------------------------------------------------
1461 subroutine read_x(kanal,*)
1463 include 'DIMENSIONS'
1464 include 'COMMON.GEO'
1465 include 'COMMON.VAR'
1466 include 'COMMON.CHAIN'
1467 include 'COMMON.IOUNITS'
1468 include 'COMMON.CONTROL'
1469 include 'COMMON.LOCAL'
1470 include 'COMMON.INTERACT'
1471 integer i,j,k,l,kanal
1472 c Read coordinates from input
1474 read(kanal,'(8f10.5)',end=10,err=10)
1475 & ((c(l,k),l=1,3),k=1,nres),
1476 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1479 c(j,2*nres)=c(j,nres)
1481 call int_from_cart1(.false.)
1484 dc(j,i)=c(j,i+1)-c(j,i)
1485 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1489 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1491 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1492 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1500 c----------------------------------------------------------------------------
1501 subroutine read_threadbase
1503 include 'DIMENSIONS'
1504 include 'COMMON.IOUNITS'
1505 include 'COMMON.GEO'
1506 include 'COMMON.VAR'
1507 include 'COMMON.INTERACT'
1508 include 'COMMON.LOCAL'
1509 include 'COMMON.NAMES'
1510 include 'COMMON.CHAIN'
1511 include 'COMMON.FFIELD'
1512 include 'COMMON.SBRIDGE'
1513 include 'COMMON.HEADER'
1514 include 'COMMON.CONTROL'
1515 include 'COMMON.DBASE'
1516 include 'COMMON.THREAD'
1517 include 'COMMON.TIME1'
1519 double precision dist
1520 C Read pattern database for threading.
1521 read (icbase,*) nseq
1523 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1524 & nres_base(2,i),nres_base(3,i)
1525 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1527 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1528 c & nres_base(2,i),nres_base(3,i)
1529 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1533 if (weidis.eq.0.0D0) weidis=0.1D0
1542 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1543 write (iout,'(a,i5)') 'nexcl: ',nexcl
1544 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1547 c------------------------------------------------------------------------------
1548 subroutine setup_var
1550 include 'DIMENSIONS'
1551 include 'COMMON.IOUNITS'
1552 include 'COMMON.GEO'
1553 include 'COMMON.VAR'
1554 include 'COMMON.INTERACT'
1555 include 'COMMON.LOCAL'
1556 include 'COMMON.NAMES'
1557 include 'COMMON.CHAIN'
1558 include 'COMMON.FFIELD'
1559 include 'COMMON.SBRIDGE'
1560 include 'COMMON.HEADER'
1561 include 'COMMON.CONTROL'
1562 include 'COMMON.DBASE'
1563 include 'COMMON.THREAD'
1564 include 'COMMON.TIME1'
1566 C Set up variable list.
1571 write (iout,*) "SETUP_VAR ialph"
1573 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1575 ialph(i,1)=nvar+nside
1579 if (indphi.gt.0) then
1581 else if (indback.gt.0) then
1586 write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1589 c----------------------------------------------------------------------------
1590 subroutine gen_dist_constr
1591 C Generate CA distance constraints.
1593 include 'DIMENSIONS'
1594 include 'COMMON.IOUNITS'
1595 include 'COMMON.GEO'
1596 include 'COMMON.VAR'
1597 include 'COMMON.INTERACT'
1598 include 'COMMON.LOCAL'
1599 include 'COMMON.NAMES'
1600 include 'COMMON.CHAIN'
1601 include 'COMMON.FFIELD'
1602 include 'COMMON.SBRIDGE'
1603 include 'COMMON.HEADER'
1604 include 'COMMON.CONTROL'
1605 include 'COMMON.DBASE'
1606 include 'COMMON.THREAD'
1607 include 'COMMON.TIME1'
1608 integer i,j,itype_pdb(maxres)
1609 common /pizda/ itype_pdb
1610 double precision dist
1612 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1613 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1614 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1616 do i=nstart_sup,nstart_sup+nsup-1
1617 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1618 cd & ' seq_pdb', restyp(itype_pdb(i))
1619 do j=i+2,nstart_sup+nsup-1
1621 ihpb(nhpb)=i+nstart_seq-nstart_sup
1622 jhpb(nhpb)=j+nstart_seq-nstart_sup
1624 dhpb(nhpb)=dist(i,j)
1627 cd write (iout,'(a)') 'Distance constraints:'
1632 cd if (ii.gt.nres) then
1637 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1638 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1639 cd & dhpb(i),forcon(i)
1643 c----------------------------------------------------------------------------
1646 include 'DIMENSIONS'
1647 include 'COMMON.MAP'
1648 include 'COMMON.IOUNITS'
1650 character*3 angid(4) /'THE','PHI','ALP','OME'/
1651 character*80 mapcard,ucase
1653 read (inp,'(a)') mapcard
1654 mapcard=ucase(mapcard)
1655 if (index(mapcard,'PHI').gt.0) then
1657 else if (index(mapcard,'THE').gt.0) then
1659 else if (index(mapcard,'ALP').gt.0) then
1661 else if (index(mapcard,'OME').gt.0) then
1664 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1665 stop 'Error - illegal variable spec in MAP card.'
1667 call readi (mapcard,'RES1',res1(imap),0)
1668 call readi (mapcard,'RES2',res2(imap),0)
1669 if (res1(imap).eq.0) then
1670 res1(imap)=res2(imap)
1671 else if (res2(imap).eq.0) then
1672 res2(imap)=res1(imap)
1674 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1676 & 'Error - illegal definition of variable group in MAP.'
1677 stop 'Error - illegal definition of variable group in MAP.'
1679 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1680 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1681 call readi(mapcard,'NSTEP',nstep(imap),0)
1682 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1684 & 'Illegal boundary and/or step size specification in MAP.'
1685 stop 'Illegal boundary and/or step size specification in MAP.'
1690 c----------------------------------------------------------------------------
1693 include 'DIMENSIONS'
1694 include 'COMMON.IOUNITS'
1695 include 'COMMON.GEO'
1696 include 'COMMON.CSA'
1697 include 'COMMON.BANK'
1698 include 'COMMON.CONTROL'
1700 character*620 mcmcard
1701 call card_concat(mcmcard)
1703 call readi(mcmcard,'NCONF',nconf,50)
1704 call readi(mcmcard,'NADD',nadd,0)
1705 call readi(mcmcard,'JSTART',jstart,1)
1706 call readi(mcmcard,'JEND',jend,1)
1707 call readi(mcmcard,'NSTMAX',nstmax,500000)
1708 call readi(mcmcard,'N0',n0,1)
1709 call readi(mcmcard,'N1',n1,6)
1710 call readi(mcmcard,'N2',n2,4)
1711 call readi(mcmcard,'N3',n3,0)
1712 call readi(mcmcard,'N4',n4,0)
1713 call readi(mcmcard,'N5',n5,0)
1714 call readi(mcmcard,'N6',n6,10)
1715 call readi(mcmcard,'N7',n7,0)
1716 call readi(mcmcard,'N8',n8,0)
1717 call readi(mcmcard,'N9',n9,0)
1718 call readi(mcmcard,'N14',n14,0)
1719 call readi(mcmcard,'N15',n15,0)
1720 call readi(mcmcard,'N16',n16,0)
1721 call readi(mcmcard,'N17',n17,0)
1722 call readi(mcmcard,'N18',n18,0)
1724 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1726 call readi(mcmcard,'NDIFF',ndiff,2)
1727 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1728 call readi(mcmcard,'IS1',is1,1)
1729 call readi(mcmcard,'IS2',is2,8)
1730 call readi(mcmcard,'NRAN0',nran0,4)
1731 call readi(mcmcard,'NRAN1',nran1,2)
1732 call readi(mcmcard,'IRR',irr,1)
1733 call readi(mcmcard,'NSEED',nseed,20)
1734 call readi(mcmcard,'NTOTAL',ntotal,10000)
1735 call reada(mcmcard,'CUT1',cut1,2.0d0)
1736 call reada(mcmcard,'CUT2',cut2,5.0d0)
1737 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1738 call readi(mcmcard,'ICMAX',icmax,3)
1739 call readi(mcmcard,'IRESTART',irestart,0)
1740 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1743 call reada(mcmcard,'DELE',dele,20.0d0)
1744 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1745 call readi(mcmcard,'IREF',iref,0)
1746 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1747 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1748 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1749 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1750 write (iout,*) "NCONF_IN",nconf_in
1753 c----------------------------------------------------------------------------
1756 include 'DIMENSIONS'
1757 include 'COMMON.MCM'
1758 include 'COMMON.MCE'
1759 include 'COMMON.IOUNITS'
1761 character*320 mcmcard
1763 call card_concat(mcmcard)
1764 call readi(mcmcard,'MAXACC',maxacc,100)
1765 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1766 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1767 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1768 call readi(mcmcard,'MAXREPM',maxrepm,200)
1769 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1770 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1771 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1772 call reada(mcmcard,'E_UP',e_up,5.0D0)
1773 call reada(mcmcard,'DELTE',delte,0.1D0)
1774 call readi(mcmcard,'NSWEEP',nsweep,5)
1775 call readi(mcmcard,'NSTEPH',nsteph,0)
1776 call readi(mcmcard,'NSTEPC',nstepc,0)
1777 call reada(mcmcard,'TMIN',tmin,298.0D0)
1778 call reada(mcmcard,'TMAX',tmax,298.0D0)
1779 call readi(mcmcard,'NWINDOW',nwindow,0)
1780 call readi(mcmcard,'PRINT_MC',print_mc,0)
1781 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1782 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1783 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1784 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1785 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1786 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1787 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1788 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1789 if (nwindow.gt.0) then
1790 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1792 winlen(i)=winend(i)-winstart(i)+1
1795 if (tmax.lt.tmin) tmax=tmin
1796 if (tmax.eq.tmin) then
1800 if (nstepc.gt.0 .and. nsteph.gt.0) then
1801 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1802 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1804 C Probabilities of different move types
1805 sumpro_type(0)=0.0D0
1806 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1807 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1808 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1809 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1810 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1811 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1812 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1814 print *,'i',i,' sumprotype',sumpro_type(i)
1815 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1816 print *,'i',i,' sumprotype',sumpro_type(i)
1820 c----------------------------------------------------------------------------
1821 subroutine read_minim
1823 include 'DIMENSIONS'
1824 include 'COMMON.MINIM'
1825 include 'COMMON.IOUNITS'
1827 character*320 minimcard
1828 call card_concat(minimcard)
1829 call readi(minimcard,'MAXMIN',maxmin,2000)
1830 call readi(minimcard,'MAXFUN',maxfun,5000)
1831 call readi(minimcard,'MINMIN',minmin,maxmin)
1832 call readi(minimcard,'MINFUN',minfun,maxmin)
1833 call reada(minimcard,'TOLF',tolf,1.0D-2)
1834 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1835 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1836 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1837 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1838 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1839 & 'Options in energy minimization:'
1840 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1841 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1842 & 'MinMin:',MinMin,' MinFun:',MinFun,
1843 & ' TolF:',TolF,' RTolF:',RTolF
1846 c----------------------------------------------------------------------------
1847 subroutine read_angles(kanal,*)
1849 include 'DIMENSIONS'
1850 include 'COMMON.GEO'
1851 include 'COMMON.VAR'
1852 include 'COMMON.CHAIN'
1853 include 'COMMON.IOUNITS'
1854 include 'COMMON.CONTROL'
1856 c Read angles from input
1858 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1859 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1860 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1861 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1864 c 9/7/01 avoid 180 deg valence angle
1865 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1867 theta(i)=deg2rad*theta(i)
1868 phi(i)=deg2rad*phi(i)
1869 alph(i)=deg2rad*alph(i)
1870 omeg(i)=deg2rad*omeg(i)
1875 c----------------------------------------------------------------------------
1876 subroutine reada(rekord,lancuch,wartosc,default)
1878 character*(*) rekord,lancuch
1879 double precision wartosc,default
1882 iread=index(rekord,lancuch)
1883 if (iread.eq.0) then
1887 iread=iread+ilen(lancuch)+1
1888 read (rekord(iread:),*,err=10,end=10) wartosc
1893 c----------------------------------------------------------------------------
1894 subroutine readi(rekord,lancuch,wartosc,default)
1896 character*(*) rekord,lancuch
1897 integer wartosc,default
1900 iread=index(rekord,lancuch)
1901 if (iread.eq.0) then
1905 iread=iread+ilen(lancuch)+1
1906 read (rekord(iread:),*,err=10,end=10) wartosc
1911 c----------------------------------------------------------------------------
1912 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1915 integer tablica(dim),default
1916 character*(*) rekord,lancuch
1923 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1924 if (iread.eq.0) return
1925 iread=iread+ilen(lancuch)+1
1926 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1929 c----------------------------------------------------------------------------
1930 subroutine multreada(rekord,lancuch,tablica,dim,default)
1933 double precision tablica(dim),default
1934 character*(*) rekord,lancuch
1941 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1942 if (iread.eq.0) return
1943 iread=iread+ilen(lancuch)+1
1944 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1947 c----------------------------------------------------------------------------
1948 subroutine openunits
1950 include 'DIMENSIONS'
1954 character*16 form,nodename
1957 include 'COMMON.SETUP'
1958 include 'COMMON.IOUNITS'
1960 include 'COMMON.CONTROL'
1961 integer lenpre,lenpot,ilen,lentmp,npos
1963 character*3 out1file_text,ucase
1968 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1969 call getenv_loc("PREFIX",prefix)
1971 call getenv_loc("POT",pot)
1972 call getenv_loc("DIRTMP",tmpdir)
1973 call getenv_loc("CURDIR",curdir)
1974 call getenv_loc("OUT1FILE",out1file_text)
1975 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1976 out1file_text=ucase(out1file_text)
1977 if (out1file_text(1:1).eq."Y") then
1980 out1file=fg_rank.gt.0
1985 if (lentmp.gt.0) then
1986 write (*,'(80(1h!))')
1987 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1988 write (*,'(80(1h!))')
1989 write (*,*)"All output files will be on node /tmp directory."
1991 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1992 if (me.eq.king) then
1993 write (*,*) "The master node is ",nodename
1994 else if (fg_rank.eq.0) then
1995 write (*,*) "I am the CG slave node ",nodename
1997 write (*,*) "I am the FG slave node ",nodename
2000 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2001 lenpre = lentmp+lenpre+1
2003 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2004 C Get the names and open the input files
2005 #if defined(WINIFL) || defined(WINPGI)
2006 open(1,file=pref_orig(:ilen(pref_orig))//
2007 & '.inp',status='old',readonly,shared)
2008 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2009 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2010 C Get parameter filenames and open the parameter files.
2011 call getenv_loc('BONDPAR',bondname)
2012 open (ibond,file=bondname,status='old',readonly,shared)
2013 call getenv_loc('THETPAR',thetname)
2014 open (ithep,file=thetname,status='old',readonly,shared)
2015 call getenv_loc('ROTPAR',rotname)
2016 open (irotam,file=rotname,status='old',readonly,shared)
2017 call getenv_loc('TORPAR',torname)
2018 open (itorp,file=torname,status='old',readonly,shared)
2019 call getenv_loc('TORDPAR',tordname)
2020 open (itordp,file=tordname,status='old',readonly,shared)
2021 call getenv_loc('FOURIER',fouriername)
2022 open (ifourier,file=fouriername,status='old',readonly,shared)
2023 call getenv_loc('ELEPAR',elename)
2024 open (ielep,file=elename,status='old',readonly,shared)
2025 call getenv_loc('SIDEPAR',sidename)
2026 open (isidep,file=sidename,status='old',readonly,shared)
2027 call getenv_loc('LIPTRANPAR',liptranname)
2028 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2029 #elif (defined CRAY) || (defined AIX)
2030 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2032 c print *,"Processor",myrank," opened file 1"
2033 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2034 c print *,"Processor",myrank," opened file 9"
2035 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2036 C Get parameter filenames and open the parameter files.
2037 call getenv_loc('BONDPAR',bondname)
2038 open (ibond,file=bondname,status='old',action='read')
2039 c print *,"Processor",myrank," opened file IBOND"
2040 call getenv_loc('THETPAR',thetname)
2041 open (ithep,file=thetname,status='old',action='read')
2043 call getenv_loc('THETPARPDB',thetname_pdb)
2044 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2046 c print *,"Processor",myrank," opened file ITHEP"
2047 call getenv_loc('ROTPAR',rotname)
2048 open (irotam,file=rotname,status='old',action='read')
2050 call getenv_loc('ROTPARPDB',rotname_pdb)
2051 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2053 c print *,"Processor",myrank," opened file IROTAM"
2054 call getenv_loc('TORPAR',torname)
2055 open (itorp,file=torname,status='old',action='read')
2056 c print *,"Processor",myrank," opened file ITORP"
2057 call getenv_loc('TORDPAR',tordname)
2058 open (itordp,file=tordname,status='old',action='read')
2059 c print *,"Processor",myrank," opened file ITORDP"
2060 call getenv_loc('SCCORPAR',sccorname)
2061 open (isccor,file=sccorname,status='old',action='read')
2062 c print *,"Processor",myrank," opened file ISCCOR"
2063 call getenv_loc('FOURIER',fouriername)
2064 open (ifourier,file=fouriername,status='old',action='read')
2065 c print *,"Processor",myrank," opened file IFOURIER"
2066 call getenv_loc('ELEPAR',elename)
2067 open (ielep,file=elename,status='old',action='read')
2068 c print *,"Processor",myrank," opened file IELEP"
2069 call getenv_loc('SIDEPAR',sidename)
2070 open (isidep,file=sidename,status='old',action='read')
2071 call getenv_loc('LIPTRANPAR',liptranname)
2072 open (iliptranpar,file=liptranname,status='old',action='read')
2073 c print *,"Processor",myrank," opened file ISIDEP"
2074 c print *,"Processor",myrank," opened parameter files"
2076 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2077 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2078 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2079 C Get parameter filenames and open the parameter files.
2080 call getenv_loc('BONDPAR',bondname)
2081 open (ibond,file=bondname,status='old')
2082 call getenv_loc('THETPAR',thetname)
2083 open (ithep,file=thetname,status='old')
2085 call getenv_loc('THETPARPDB',thetname_pdb)
2086 open (ithep_pdb,file=thetname_pdb,status='old')
2088 call getenv_loc('ROTPAR',rotname)
2089 open (irotam,file=rotname,status='old')
2091 call getenv_loc('ROTPARPDB',rotname_pdb)
2092 open (irotam_pdb,file=rotname_pdb,status='old')
2094 call getenv_loc('TORPAR',torname)
2095 open (itorp,file=torname,status='old')
2096 call getenv_loc('TORDPAR',tordname)
2097 open (itordp,file=tordname,status='old')
2098 call getenv_loc('SCCORPAR',sccorname)
2099 open (isccor,file=sccorname,status='old')
2100 call getenv_loc('FOURIER',fouriername)
2101 open (ifourier,file=fouriername,status='old')
2102 call getenv_loc('ELEPAR',elename)
2103 open (ielep,file=elename,status='old')
2104 call getenv_loc('SIDEPAR',sidename)
2105 open (isidep,file=sidename,status='old')
2106 call getenv_loc('LIPTRANPAR',liptranname)
2107 open (iliptranpar,file=liptranname,status='old')
2109 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2111 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2112 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2113 C Get parameter filenames and open the parameter files.
2114 call getenv_loc('BONDPAR',bondname)
2115 open (ibond,file=bondname,status='old',readonly)
2116 call getenv_loc('THETPAR',thetname)
2117 open (ithep,file=thetname,status='old',readonly)
2118 call getenv_loc('ROTPAR',rotname)
2119 open (irotam,file=rotname,status='old',readonly)
2120 call getenv_loc('TORPAR',torname)
2121 open (itorp,file=torname,status='old',readonly)
2122 call getenv_loc('TORDPAR',tordname)
2123 open (itordp,file=tordname,status='old',readonly)
2124 call getenv_loc('SCCORPAR',sccorname)
2125 open (isccor,file=sccorname,status='old',readonly)
2127 call getenv_loc('THETPARPDB',thetname_pdb)
2128 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2130 call getenv_loc('FOURIER',fouriername)
2131 open (ifourier,file=fouriername,status='old',readonly)
2132 call getenv_loc('ELEPAR',elename)
2133 open (ielep,file=elename,status='old',readonly)
2134 call getenv_loc('SIDEPAR',sidename)
2135 open (isidep,file=sidename,status='old',readonly)
2136 call getenv_loc('LIPTRANPAR',liptranname)
2137 open (iliptranpar,file=liptranname,status='old',action='read')
2139 call getenv_loc('ROTPARPDB',rotname_pdb)
2140 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2144 call getenv_loc('TUBEPAR',tubename)
2145 #if defined(WINIFL) || defined(WINPGI)
2146 open (itube,file=tubename,status='old',readonly,shared)
2147 #elif (defined CRAY) || (defined AIX)
2148 open (itube,file=tubename,status='old',action='read')
2150 open (itube,file=tubename,status='old')
2152 open (itube,file=tubename,status='old',readonly)
2157 C 8/9/01 In the newest version SCp interaction constants are read from a file
2158 C Use -DOLDSCP to use hard-coded constants instead.
2160 call getenv_loc('SCPPAR',scpname)
2161 #if defined(WINIFL) || defined(WINPGI)
2162 open (iscpp,file=scpname,status='old',readonly,shared)
2163 #elif (defined CRAY) || (defined AIX)
2164 open (iscpp,file=scpname,status='old',action='read')
2166 open (iscpp,file=scpname,status='old')
2168 open (iscpp,file=scpname,status='old',readonly)
2171 call getenv_loc('PATTERN',patname)
2172 #if defined(WINIFL) || defined(WINPGI)
2173 open (icbase,file=patname,status='old',readonly,shared)
2174 #elif (defined CRAY) || (defined AIX)
2175 open (icbase,file=patname,status='old',action='read')
2177 open (icbase,file=patname,status='old')
2179 open (icbase,file=patname,status='old',readonly)
2182 C Open output file only for CG processes
2183 c print *,"Processor",myrank," fg_rank",fg_rank
2184 if (fg_rank.eq.0) then
2186 if (nodes.eq.1) then
2189 npos = dlog10(dfloat(nodes-1))+1
2191 if (npos.lt.3) npos=3
2192 write (liczba,'(i1)') npos
2193 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2195 write (liczba,form) me
2196 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2197 & liczba(:ilen(liczba))
2198 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2200 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2202 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2203 & liczba(:ilen(liczba))//'.mol2'
2204 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2205 & liczba(:ilen(liczba))//'.stat'
2207 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2208 & //liczba(:ilen(liczba))//'.stat')
2209 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2212 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2213 & liczba(:ilen(liczba))//'.const'
2218 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2219 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2220 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2221 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2222 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2224 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2226 rest2name=prefix(:ilen(prefix))//'.rst'
2228 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2231 #if defined(AIX) || defined(PGI) || defined(CRAY)
2232 if (me.eq.king .or. .not. out1file) then
2233 open(iout,file=outname,status='unknown')
2235 open(iout,file="/dev/null",status="unknown")
2239 if (fg_rank.gt.0) then
2240 write (liczba,'(i3.3)') myrank/nfgtasks
2241 write (ll,'(bz,i3.3)') fg_rank
2242 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2248 open(igeom,file=intname,status='unknown',position='append')
2249 open(ipdb,file=pdbname,status='unknown')
2250 open(imol2,file=mol2name,status='unknown')
2251 open(istat,file=statname,status='unknown',position='append')
2253 c1out open(iout,file=outname,status='unknown')
2256 if (me.eq.king .or. .not.out1file)
2257 & open(iout,file=outname,status='unknown')
2259 if (fg_rank.gt.0) then
2260 write (liczba,'(i3.3)') myrank/nfgtasks
2261 write (ll,'(bz,i3.3)') fg_rank
2262 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2267 open(igeom,file=intname,status='unknown',access='append')
2268 open(ipdb,file=pdbname,status='unknown')
2269 open(imol2,file=mol2name,status='unknown')
2270 open(istat,file=statname,status='unknown',access='append')
2272 c1out open(iout,file=outname,status='unknown')
2275 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2276 csa_seed=prefix(:lenpre)//'.CSA.seed'
2277 csa_history=prefix(:lenpre)//'.CSA.history'
2278 csa_bank=prefix(:lenpre)//'.CSA.bank'
2279 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2280 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2281 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2282 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2283 csa_int=prefix(:lenpre)//'.int'
2284 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2285 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2286 csa_in=prefix(:lenpre)//'.CSA.in'
2287 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2290 write (iout,'(80(1h-))')
2291 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2292 write (iout,'(80(1h-))')
2293 write (iout,*) "Input file : ",
2294 & pref_orig(:ilen(pref_orig))//'.inp'
2295 write (iout,*) "Output file : ",
2296 & outname(:ilen(outname))
2298 write (iout,*) "Sidechain potential file : ",
2299 & sidename(:ilen(sidename))
2301 write (iout,*) "SCp potential file : ",
2302 & scpname(:ilen(scpname))
2304 write (iout,*) "Electrostatic potential file : ",
2305 & elename(:ilen(elename))
2306 write (iout,*) "Cumulant coefficient file : ",
2307 & fouriername(:ilen(fouriername))
2308 write (iout,*) "Torsional parameter file : ",
2309 & torname(:ilen(torname))
2310 write (iout,*) "Double torsional parameter file : ",
2311 & tordname(:ilen(tordname))
2312 write (iout,*) "SCCOR parameter file : ",
2313 & sccorname(:ilen(sccorname))
2314 write (iout,*) "Bond & inertia constant file : ",
2315 & bondname(:ilen(bondname))
2316 write (iout,*) "Bending-torsion parameter file : ",
2317 & thetname(:ilen(thetname))
2318 write (iout,*) "Rotamer parameter file : ",
2319 & rotname(:ilen(rotname))
2320 write (iout,*) "Thetpdb parameter file : ",
2321 & thetname_pdb(:ilen(thetname_pdb))
2322 write (iout,*) "Threading database : ",
2323 & patname(:ilen(patname))
2325 &write (iout,*)" DIRTMP : ",
2327 write (iout,'(80(1h-))')
2331 c----------------------------------------------------------------------------
2332 subroutine card_concat(card)
2333 implicit real*8 (a-h,o-z)
2334 include 'DIMENSIONS'
2335 include 'COMMON.IOUNITS'
2337 character*80 karta,ucase
2339 read (inp,'(a)') karta
2342 do while (karta(80:80).eq.'&')
2343 card=card(:ilen(card)+1)//karta(:79)
2344 read (inp,'(a)') karta
2347 card=card(:ilen(card)+1)//karta
2350 c------------------------------------------------------------------------------
2353 include 'DIMENSIONS'
2354 include 'COMMON.CHAIN'
2355 include 'COMMON.IOUNITS'
2356 include 'COMMON.CONTROL'
2358 include 'COMMON.QRESTR'
2360 include 'COMMON.LAGRANGE.5diag'
2362 include 'COMMON.LAGRANGE'
2365 open(irest2,file=rest2name,status='unknown')
2366 read(irest2,*) totT,EK,potE,totE,t_bath
2369 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2372 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2375 read (irest2,*) iset
2380 c------------------------------------------------------------------------------
2381 subroutine read_fragments
2383 include 'DIMENSIONS'
2387 include 'COMMON.SETUP'
2388 include 'COMMON.CHAIN'
2389 include 'COMMON.IOUNITS'
2391 include 'COMMON.QRESTR'
2392 include 'COMMON.CONTROL'
2394 read(inp,*) nset,nfrag,npair,nfrag_back
2395 loc_qlike=(nfrag_back.lt.0)
2396 nfrag_back=iabs(nfrag_back)
2397 if(me.eq.king.or..not.out1file)
2398 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2399 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2401 read(inp,*) mset(iset)
2403 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2405 if(me.eq.king.or..not.out1file)
2406 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2407 & ifrag(2,i,iset), qinfrag(i,iset)
2410 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2412 if(me.eq.king.or..not.out1file)
2413 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2414 & ipair(2,i,iset), qinpair(i,iset)
2418 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2419 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2420 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2421 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2422 if(me.eq.king.or..not.out1file)
2423 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2424 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2425 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2426 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2430 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2431 & wfrag_back(3,i,iset),
2432 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2433 if(me.eq.king.or..not.out1file)
2434 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2435 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2441 C---------------------------------------------------------------------------
2442 subroutine read_afminp
2444 include 'DIMENSIONS'
2448 include 'COMMON.SETUP'
2449 include 'COMMON.CONTROL'
2450 include 'COMMON.CHAIN'
2451 include 'COMMON.IOUNITS'
2452 include 'COMMON.SBRIDGE'
2453 character*320 afmcard
2455 c print *, "wchodze"
2456 call card_concat(afmcard)
2457 call readi(afmcard,"BEG",afmbeg,0)
2458 call readi(afmcard,"END",afmend,0)
2459 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2460 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2461 print *,'FORCE=' ,forceAFMconst
2462 CCCC NOW PROPERTIES FOR AFM
2465 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2467 distafminit=dsqrt(distafminit)
2468 c print *,'initdist',distafminit
2471 c-------------------------------------------------------------------------------
2472 subroutine read_saxs_constr
2474 include 'DIMENSIONS'
2478 include 'COMMON.SETUP'
2479 include 'COMMON.CONTROL'
2480 include 'COMMON.SAXS'
2481 include 'COMMON.CHAIN'
2482 include 'COMMON.IOUNITS'
2483 include 'COMMON.SBRIDGE'
2484 double precision cm(3),cnorm
2487 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2489 if (saxs_mode.eq.0) then
2490 c SAXS distance distribution
2492 read(inp,*) distsaxs(i),Psaxs(i)
2496 Cnorm = Cnorm + Psaxs(i)
2498 write (iout,*) "Cnorm",Cnorm
2500 Psaxs(i)=Psaxs(i)/Cnorm
2502 write (iout,*) "Normalized distance distribution from SAXS"
2504 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2508 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2510 write (iout,*) "Wsaxs0",Wsaxs0
2514 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2521 cm(j)=cm(j)+Csaxs(j,i)
2529 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2532 write (iout,*) "SAXS sphere coordinates"
2534 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2540 c-------------------------------------------------------------------------------
2541 subroutine read_dist_constr
2543 include 'DIMENSIONS'
2547 include 'COMMON.SETUP'
2548 include 'COMMON.CONTROL'
2549 include 'COMMON.CHAIN'
2550 include 'COMMON.IOUNITS'
2551 include 'COMMON.SBRIDGE'
2552 include 'COMMON.INTERACT'
2553 integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2554 integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2555 double precision wfrag_(100),wpair_(1000)
2556 double precision ddjk,dist,dist_cut,fordepthmax
2557 character*5000 controlcard
2558 logical normalize,next
2560 double precision scal_bfac
2561 double precision xlink(4,0:4) /
2563 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2564 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2565 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2566 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2567 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2568 c print *, "WCHODZE"
2569 write (iout,*) "Calling read_dist_constr"
2570 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2572 restr_on_coord=.false.
2581 call card_concat(controlcard)
2582 next = index(controlcard,"NEXT").gt.0
2583 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2584 write (iout,*) "restr_type",restr_type
2585 call readi(controlcard,"NFRAG",nfrag_,0)
2586 call readi(controlcard,"NFRAG",nfrag_,0)
2587 call readi(controlcard,"NPAIR",npair_,0)
2588 call readi(controlcard,"NDIST",ndist_,0)
2589 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2590 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2591 if (restr_type.eq.10)
2592 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2593 if (restr_type.eq.12)
2594 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2595 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2596 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2597 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2598 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2599 normalize = index(controlcard,"NORMALIZE").gt.0
2600 write (iout,*) "WBOLTZD",wboltzd
2601 write (iout,*) "SCAL_PEAK",scal_peak
2602 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2603 write (iout,*) "IFRAG"
2605 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2607 write (iout,*) "IPAIR"
2609 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2611 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2613 & "Distance restraints as generated from reference structure"
2615 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2616 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2617 & ifrag_(2,i)=nstart_sup+nsup-1
2618 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2620 if (wfrag_(i).eq.0.0d0) cycle
2621 do j=ifrag_(1,i),ifrag_(2,i)-1
2622 do k=j+1,ifrag_(2,i)
2623 c write (iout,*) "j",j," k",k
2625 if (restr_type.eq.1) then
2631 forcon(nhpb)=wfrag_(i)
2632 else if (constr_dist.eq.2) then
2633 if (ddjk.le.dist_cut) then
2639 forcon(nhpb)=wfrag_(i)
2641 else if (restr_type.eq.3) then
2647 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2650 if (.not.out1file .or. me.eq.king)
2651 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2652 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2654 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2655 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2661 if (wpair_(i).eq.0.0d0) cycle
2669 do j=ifrag_(1,ii),ifrag_(2,ii)
2670 do k=ifrag_(1,jj),ifrag_(2,jj)
2672 if (restr_type.eq.1) then
2678 forcon(nhpb)=wpair_(i)
2679 else if (constr_dist.eq.2) then
2680 if (ddjk.le.dist_cut) then
2686 forcon(nhpb)=wpair_(i)
2688 else if (restr_type.eq.3) then
2694 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2697 if (.not.out1file .or. me.eq.king)
2698 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2699 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2701 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2702 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2709 write (iout,*) "Distance restraints as read from input"
2711 if (restr_type.eq.12) then
2712 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2713 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2714 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2715 & fordepth_peak(nhpb_peak+1),npeak
2716 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2717 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2718 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2719 c & fordepth_peak(nhpb_peak+1),npeak
2720 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2721 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2722 nhpb_peak=nhpb_peak+1
2723 irestr_type_peak(nhpb_peak)=12
2724 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2727 if (.not.out1file .or. me.eq.king)
2728 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2729 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2730 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2731 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2732 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2734 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2735 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2736 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2737 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2738 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2740 if (ibecarb_peak(nhpb_peak).eq.3) then
2741 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2742 else if (ibecarb_peak(nhpb_peak).eq.2) then
2743 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2744 else if (ibecarb_peak(nhpb_peak).eq.1) then
2745 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2746 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2748 else if (restr_type.eq.11) then
2749 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2750 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2751 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2752 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2754 irestr_type(nhpb)=11
2756 if (.not.out1file .or. me.eq.king)
2757 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2758 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2759 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2761 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2762 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2763 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2765 c if (ibecarb(nhpb).gt.0) then
2766 c ihpb(nhpb)=ihpb(nhpb)+nres
2767 c jhpb(nhpb)=jhpb(nhpb)+nres
2769 if (ibecarb(nhpb).eq.3) then
2770 ihpb(nhpb)=ihpb(nhpb)+nres
2771 else if (ibecarb(nhpb).eq.2) then
2772 ihpb(nhpb)=ihpb(nhpb)+nres
2773 else if (ibecarb(nhpb).eq.1) 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).eq.3) then
2812 jhpb(nhpb)=jhpb(nhpb)+nres
2813 else if (ibecarb(nhpb).eq.2) then
2814 ihpb(nhpb)=ihpb(nhpb)+nres
2815 else if (ibecarb(nhpb).eq.1) then
2816 ihpb(nhpb)=ihpb(nhpb)+nres
2817 jhpb(nhpb)=jhpb(nhpb)+nres
2820 if (.not.out1file .or. me.eq.king)
2821 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2822 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2823 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2826 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2827 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2828 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2833 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2834 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2835 if (forcon(nhpb+1).gt.0.0d0) then
2837 if (dhpb1(nhpb).eq.0.0d0) then
2842 if (ibecarb(nhpb).eq.3) then
2843 jhpb(nhpb)=jhpb(nhpb)+nres
2844 else if (ibecarb(nhpb).eq.2) then
2845 ihpb(nhpb)=ihpb(nhpb)+nres
2846 else if (ibecarb(nhpb).eq.1) then
2847 ihpb(nhpb)=ihpb(nhpb)+nres
2848 jhpb(nhpb)=jhpb(nhpb)+nres
2850 if (dhpb(nhpb).eq.0.0d0)
2851 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2854 if (.not.out1file .or. me.eq.king)
2855 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2856 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2858 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2859 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2862 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2863 C if (forcon(nhpb+1).gt.0.0d0) then
2865 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2868 if (restr_type.eq.4) then
2869 write (iout,*) "The BFAC array"
2871 write (iout,'(i5,f10.5)') i,bfac(i)
2874 if (itype(i).eq.ntyp1) cycle
2876 if (itype(j).eq.ntyp1) cycle
2877 if (itype(i).eq.10) then
2882 if (itype(j).eq.10) then
2892 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2895 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2896 ihpb(nhpb)=i+nres*ii
2897 jhpb(nhpb)=j+nres*jj
2898 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2900 if (.not.out1file .or. me.eq.king) then
2901 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2902 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2903 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2907 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2908 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2909 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2919 if (restr_type.eq.5) then
2920 restr_on_coord=.true.
2922 if (itype(i).eq.ntyp1) cycle
2923 bfac(i)=(scal_bfac/bfac(i))**2
2932 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
2933 & fordepthmax=fordepth(i)
2936 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
2941 c-------------------------------------------------------------------------------
2942 subroutine read_constr_homology
2944 include 'DIMENSIONS'
2948 include 'COMMON.SETUP'
2949 include 'COMMON.CONTROL'
2950 include 'COMMON.HOMOLOGY'
2951 include 'COMMON.CHAIN'
2952 include 'COMMON.IOUNITS'
2954 include 'COMMON.QRESTR'
2955 include 'COMMON.GEO'
2956 include 'COMMON.INTERACT'
2957 include 'COMMON.NAMES'
2959 c For new homol impl
2961 include 'COMMON.VAR'
2964 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2966 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2967 c & sigma_odl_temp(maxres,maxres,max_template)
2969 character*24 model_ki_dist, model_ki_angle
2970 character*500 controlcard
2971 integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
2972 & ik,iistart,iishift
2977 c FP - Nov. 2014 Temporary specifications for new vars
2979 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2981 double precision, dimension (max_template,maxres) :: rescore
2982 double precision, dimension (max_template,maxres) :: rescore2
2983 double precision, dimension (max_template,maxres) :: rescore3
2984 double precision distal
2985 character*24 pdbfile,tpl_k_rescore
2986 c -----------------------------------------------------------------
2987 c Reading multiple PDB ref structures and calculation of retraints
2988 c not using pre-computed ones stored in files model_ki_{dist,angle}
2990 c -----------------------------------------------------------------
2993 c Alternative: reading from input
2994 call card_concat(controlcard)
2995 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2996 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2997 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2998 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2999 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3000 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3001 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3002 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3003 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3004 if(.not.read2sigma.and.start_from_model) then
3005 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3006 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3007 start_from_model=.false.
3009 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3010 & write(iout,*) 'START_FROM_MODELS is ON'
3011 if(start_from_model .and. rest) then
3012 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3013 write(iout,*) 'START_FROM_MODELS is OFF'
3014 write(iout,*) 'remove restart keyword from input'
3017 if (homol_nset.gt.1)then
3018 call card_concat(controlcard)
3019 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3020 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3021 write(iout,*) "iset homology_weight "
3023 write(iout,*) i,waga_homology(i)
3026 iset=mod(kolor,homol_nset)+1
3029 waga_homology(1)=1.0
3032 cd write (iout,*) "nnt",nnt," nct",nct
3039 c write(iout,*) 'nnt=',nnt,'nct=',nct
3042 do k=1,constr_homology
3055 if (read_homol_frag) then
3056 call read_klapaucjusz
3059 do k=1,constr_homology
3061 read(inp,'(a)') pdbfile
3062 if(me.eq.king .or. .not. out1file)
3063 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3064 & pdbfile(:ilen(pdbfile))
3065 open(ipdbin,file=pdbfile,status='old',err=33)
3067 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3068 & pdbfile(:ilen(pdbfile))
3071 c print *,'Begin reading pdb data'
3073 c Files containing res sim or local scores (former containing sigmas)
3076 write(kic2,'(bz,i2.2)') k
3078 tpl_k_rescore="template"//kic2//".sco"
3081 if (read2sigma) then
3082 call readpdb_template(k)
3087 c Distance restraints
3090 C Copy the coordinates from reference coordinates (?)
3094 c write (iout,*) "c(",j,i,") =",c(j,i)
3098 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3100 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3101 open (ientin,file=tpl_k_rescore,status='old')
3102 if (nnt.gt.1) rescore(k,1)=0.0d0
3103 do irec=nnt,nct ! loop for reading res sim
3104 if (read2sigma) then
3105 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3106 & rescore3_tmp,idomain_tmp
3108 idomain(k,i_tmp)=idomain_tmp
3109 rescore(k,i_tmp)=rescore_tmp
3110 rescore2(k,i_tmp)=rescore2_tmp
3111 rescore3(k,i_tmp)=rescore3_tmp
3112 if (.not. out1file .or. me.eq.king)
3113 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3114 & i_tmp,rescore2_tmp,rescore_tmp,
3115 & rescore3_tmp,idomain_tmp
3118 read (ientin,*,end=1401) rescore_tmp
3120 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3121 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3122 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3127 if (waga_dist.ne.0.0d0) then
3135 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3136 c write (iout,*) k,i,j,distal,dist2_cut
3138 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3139 & .and. distal.le.dist2_cut ) then
3145 c write (iout,*) "k",k
3146 c write (iout,*) "i",i," j",j," constr_homology",
3151 if (read2sigma) then
3154 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3156 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3157 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3158 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3160 if (odl(k,ii).le.dist_cut) then
3161 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3164 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3165 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3167 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3168 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3172 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3175 l_homo(k,ii)=.false.
3182 c Theta, dihedral and SC retraints
3184 if (waga_angle.gt.0.0d0) then
3185 c open (ientin,file=tpl_k_sigma_dih,status='old')
3186 c do irec=1,maxres-3 ! loop for reading sigma_dih
3187 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3188 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3189 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3190 c & sigma_dih(k,i+nnt-1)
3195 if (idomain(k,i).eq.0) then
3199 dih(k,i)=phiref(i) ! right?
3200 c read (ientin,*) sigma_dih(k,i) ! original variant
3201 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3202 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3203 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3204 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3205 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3207 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3208 & rescore(k,i-2)+rescore(k,i-3))/4.0
3209 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3210 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3211 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3212 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3213 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3214 c Instead of res sim other local measure of b/b str reliability possible
3215 if (sigma_dih(k,i).ne.0)
3216 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3217 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3222 if (waga_theta.gt.0.0d0) then
3223 c open (ientin,file=tpl_k_sigma_theta,status='old')
3224 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3225 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3226 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3227 c & sigma_theta(k,i+nnt-1)
3232 do i = nnt+2,nct ! right? without parallel.
3233 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3234 c do i=ithet_start,ithet_end ! with FG parallel.
3235 if (idomain(k,i).eq.0) then
3236 sigma_theta(k,i)=0.0
3239 thetatpl(k,i)=thetaref(i)
3240 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3241 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3242 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3243 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3244 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3245 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3246 & rescore(k,i-2))/3.0
3247 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3248 if (sigma_theta(k,i).ne.0)
3249 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3251 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3252 c rescore(k,i-2) ! right expression ?
3253 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3257 if (waga_d.gt.0.0d0) then
3258 c open (ientin,file=tpl_k_sigma_d,status='old')
3259 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3260 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3261 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3262 c & sigma_d(k,i+nnt-1)
3266 do i = nnt,nct ! right? without parallel.
3267 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3268 c do i=loc_start,loc_end ! with FG parallel.
3269 if (itype(i).eq.10) cycle
3270 if (idomain(k,i).eq.0 ) then
3277 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3278 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3279 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3280 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3281 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3282 if (sigma_d(k,i).ne.0)
3283 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3285 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3286 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3287 c read (ientin,*) sigma_d(k,i) ! 1st variant
3292 c remove distance restraints not used in any model from the list
3293 c shift data in all arrays
3295 if (waga_dist.ne.0.0d0) then
3301 if (ii_in_use(ii).eq.0.and.liiflag) then
3305 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3306 & .not.liiflag.and.ii.eq.lim_odl) then
3307 if (ii.eq.lim_odl) then
3308 iishift=ii-iistart+1
3313 do ki=iistart,lim_odl-iishift
3314 ires_homo(ki)=ires_homo(ki+iishift)
3315 jres_homo(ki)=jres_homo(ki+iishift)
3316 ii_in_use(ki)=ii_in_use(ki+iishift)
3317 do k=1,constr_homology
3318 odl(k,ki)=odl(k,ki+iishift)
3319 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3320 l_homo(k,ki)=l_homo(k,ki+iishift)
3324 lim_odl=lim_odl-iishift
3330 endif ! .not. klapaucjusz
3332 if (constr_homology.gt.0) call homology_partition
3333 if (constr_homology.gt.0) call init_int_table
3334 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3335 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3339 if (.not.out_template_restr) return
3340 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3341 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3342 write (iout,*) "Distance restraints from templates"
3344 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3345 & ii,ires_homo(ii),jres_homo(ii),
3346 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3347 & ki=1,constr_homology)
3349 write (iout,*) "Dihedral angle restraints from templates"
3351 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3352 & (rad2deg*dih(ki,i),
3353 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3355 write (iout,*) "Virtual-bond angle restraints from templates"
3357 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3358 & (rad2deg*thetatpl(ki,i),
3359 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3361 write (iout,*) "SC restraints from templates"
3363 write(iout,'(i5,100(4f8.2,4x))') i,
3364 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3365 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3368 c -----------------------------------------------------------------
3371 c----------------------------------------------------------------------
3373 subroutine flush(iu)
3378 subroutine flush(iu)
3383 c------------------------------------------------------------------------------
3384 subroutine copy_to_tmp(source)
3386 include "DIMENSIONS"
3387 include "COMMON.IOUNITS"
3388 character*(*) source
3389 character* 256 tmpfile
3393 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3394 inquire(file=tmpfile,exist=ex)
3396 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3397 & " to temporary directory..."
3398 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3399 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3403 c------------------------------------------------------------------------------
3404 subroutine move_from_tmp(source)
3406 include "DIMENSIONS"
3407 include "COMMON.IOUNITS"
3408 character*(*) source
3411 write (*,*) "Moving ",source(:ilen(source)),
3412 & " from temporary directory to working directory"
3413 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3414 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3417 c------------------------------------------------------------------------------
3418 subroutine random_init(seed)
3420 C Initialize random number generator
3423 include 'DIMENSIONS'
3426 logical OKRandom, prng_restart
3428 integer iseed_array(4)
3429 integer error_msg,ierr
3431 include 'COMMON.IOUNITS'
3432 include 'COMMON.TIME1'
3433 include 'COMMON.THREAD'
3434 include 'COMMON.SBRIDGE'
3435 include 'COMMON.CONTROL'
3436 include 'COMMON.MCM'
3437 include 'COMMON.MAP'
3438 include 'COMMON.HEADER'
3439 include 'COMMON.CSA'
3440 include 'COMMON.CHAIN'
3441 include 'COMMON.MUCA'
3443 include 'COMMON.FFIELD'
3444 include 'COMMON.SETUP'
3446 double precision seed,ran_number
3447 iseed=-dint(dabs(seed))
3448 if (iseed.eq.0) then
3449 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3450 & 'Random seed undefined. The program will stop.'
3451 write (*,'(/80(1h*)/20x,a/80(1h*))')
3452 & 'Random seed undefined. The program will stop.'
3454 call mpi_finalize(mpi_comm_world,ierr)
3456 stop 'Bad random seed.'
3459 if (fg_rank.eq.0) then
3463 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3464 OKRandom = prng_restart(me,iseed)
3467 tmp=65536.0d0**(4-i)
3468 iseed_array(i) = dint(seed/tmp)
3469 seed=seed-iseed_array(i)*tmp
3472 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3473 & (iseed_array(i),i=1,4)
3474 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3475 & (iseed_array(i),i=1,4)
3476 OKRandom = prng_restart(me,iseed_array)
3479 c r1 = prng_next(me)
3480 r1=ran_number(0.0D0,1.0D0)
3482 & write (iout,*) 'ran_num',r1
3483 if (r1.lt.0.0d0) OKRandom=.false.
3485 if (.not.OKRandom) then
3486 write (iout,*) 'PRNG IS NOT WORKING!!!'
3487 print *,'PRNG IS NOT WORKING!!!'
3490 call mpi_abort(mpi_comm_world,error_msg,ierr)
3493 write (iout,*) 'too many processors for parallel prng'
3494 write (*,*) 'too many processors for parallel prng'
3502 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3506 c----------------------------------------------------------------------
3507 subroutine read_klapaucjusz
3509 include 'DIMENSIONS'
3513 include 'COMMON.SETUP'
3514 include 'COMMON.CONTROL'
3515 include 'COMMON.HOMOLOGY'
3516 include 'COMMON.CHAIN'
3517 include 'COMMON.IOUNITS'
3519 include 'COMMON.GEO'
3520 include 'COMMON.INTERACT'
3521 include 'COMMON.NAMES'
3522 character*256 fragfile
3523 integer ninclust(maxclust),inclust(max_template,maxclust),
3524 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
3527 character*24 model_ki_dist, model_ki_angle
3528 character*500 controlcard
3529 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3530 & ik,ll,ii,kk,iistart,iishift,lim_xx
3531 double precision distal
3532 logical lprn /.true./
3538 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3539 double precision, dimension (max_template,maxres) :: rescore
3540 double precision, dimension (max_template,maxres) :: rescore2
3541 character*24 pdbfile,tpl_k_rescore
3544 c For new homol impl
3546 include 'COMMON.VAR'
3548 call getenv("FRAGFILE",fragfile)
3549 open(ientin,file=fragfile,status="old",err=10)
3550 read(ientin,*) constr_homology,nclust
3556 do k=1,constr_homology
3557 read(ientin,'(a)') pdbfile
3558 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3559 & pdbfile(:ilen(pdbfile))
3560 open(ipdbin,file=pdbfile,status='old',err=33)
3562 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3563 & pdbfile(:ilen(pdbfile))
3567 call readpdb_template(k)
3575 read(ientin,*) ninclust(i),nresclust(i)
3576 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3577 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3580 c Loop over clusters
3583 do ll = 1,ninclust(l)
3591 idomain(k,iresclust(i,l)+1) = 1
3593 idomain(k,iresclust(i,l)) = 1
3597 c Distance restraints
3600 C Copy the coordinates from reference coordinates (?)
3604 c write (iout,*) "c(",j,i,") =",c(j,i)
3607 call int_from_cart(.true.,.false.)
3608 call sc_loc_geom(.true.)
3610 thetaref(i)=theta(i)
3613 if (waga_dist.ne.0.0d0) then
3621 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3622 c write (iout,*) k,i,j,distal,dist2_cut
3624 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3625 & .and. distal.le.dist2_cut ) then
3631 c write (iout,*) "k",k
3632 c write (iout,*) "i",i," j",j," constr_homology",
3637 if (read2sigma) then
3640 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3642 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3643 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3644 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3646 if (odl(k,ii).le.dist_cut) then
3647 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3650 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3651 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3653 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3654 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3658 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3661 c l_homo(k,ii)=.false.
3668 c Theta, dihedral and SC retraints
3670 if (waga_angle.gt.0.0d0) then
3672 if (idomain(k,i).eq.0) then
3673 c sigma_dih(k,i)=0.0
3677 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3678 & rescore(k,i-2)+rescore(k,i-3))/4.0
3679 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3680 c & " sigma_dihed",sigma_dih(k,i)
3681 if (sigma_dih(k,i).ne.0)
3682 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3687 if (waga_theta.gt.0.0d0) then
3689 if (idomain(k,i).eq.0) then
3690 c sigma_theta(k,i)=0.0
3693 thetatpl(k,i)=thetaref(i)
3694 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3695 & rescore(k,i-2))/3.0
3696 if (sigma_theta(k,i).ne.0)
3697 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3701 if (waga_d.gt.0.0d0) then
3703 if (itype(i).eq.10) cycle
3704 if (idomain(k,i).eq.0 ) then
3711 sigma_d(k,i)=rescore(k,i)
3712 if (sigma_d(k,i).ne.0)
3713 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3714 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3720 c remove distance restraints not used in any model from the list
3721 c shift data in all arrays
3723 if (waga_dist.ne.0.0d0) then
3729 if (ii_in_use(ii).eq.0.and.liiflag) then
3733 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3734 & .not.liiflag.and.ii.eq.lim_odl) then
3735 if (ii.eq.lim_odl) then
3736 iishift=ii-iistart+1
3741 do ki=iistart,lim_odl-iishift
3742 ires_homo(ki)=ires_homo(ki+iishift)
3743 jres_homo(ki)=jres_homo(ki+iishift)
3744 ii_in_use(ki)=ii_in_use(ki+iishift)
3745 do k=1,constr_homology
3746 odl(k,ki)=odl(k,ki+iishift)
3747 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3748 l_homo(k,ki)=l_homo(k,ki+iishift)
3752 lim_odl=lim_odl-iishift
3759 10 stop "Error in fragment file"