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 readi(controlcard,"INTER_LIST_UPDATE",imatupdate,100)
175 call reada(controlcard,"T_BATH",t_bath,300.0d0)
176 minim=(index(controlcard,'MINIMIZE').gt.0)
177 dccart=(index(controlcard,'CART').gt.0)
178 overlapsc=(index(controlcard,'OVERLAP').gt.0)
179 overlapsc=.not.overlapsc
180 searchsc=(index(controlcard,'SEARCHSC').gt.0 .and.
181 & index(controlcard,'NOSEARCHSC').eq.0)
182 sideadd=(index(controlcard,'SIDEADD').gt.0)
183 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
184 mremd_dec=(index(controlcard,'MREMD_DEC').gt.0)
185 outpdb=(index(controlcard,'PDBOUT').gt.0)
186 outmol2=(index(controlcard,'MOL2OUT').gt.0)
187 pdbref=(index(controlcard,'PDBREF').gt.0)
188 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
189 indpdb=index(controlcard,'PDBSTART')
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 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
299 extconf=(index(controlcard,'EXTCONF').gt.0)
300 if (start_from_model) then
304 i2ndstr=index(controlcard,'USE_SEC_PRED')
305 gradout=index(controlcard,'GRADOUT').gt.0
306 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
307 C DISTCHAINMAX become obsolete for periodic boundry condition
308 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
309 C Reading the dimensions of box in x,y,z coordinates
310 call reada(controlcard,'BOXX',boxxsize,100.0d0)
311 call reada(controlcard,'BOXY',boxysize,100.0d0)
312 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
313 write(iout,*) "Periodic box dimensions",boxxsize,boxysize,boxzsize
314 c Cutoff range for interactions
315 call reada(controlcard,"R_CUT_INT",r_cut_int,25.0d0)
316 call reada(controlcard,"R_CUT_RESPA",r_cut_respa,2.0d0)
317 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
318 write (iout,*) "Cutoff on interactions",r_cut_int
320 & "Cutoff in switching short and long range interactions in RESPA",
322 write (iout,*) "lambda in switch function",rlamb
323 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
324 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
325 if (lipthick.gt.0.0d0) then
326 bordliptop=(boxzsize+lipthick)/2.0
327 bordlipbot=bordliptop-lipthick
329 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
330 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
331 buflipbot=bordlipbot+lipbufthick
332 bufliptop=bordliptop-lipbufthick
333 if ((lipbufthick*2.0d0).gt.lipthick)
334 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
336 write(iout,*) "bordliptop=",bordliptop
337 write(iout,*) "bordlipbot=",bordlipbot
338 write(iout,*) "bufliptop=",bufliptop
339 write(iout,*) "buflipbot=",buflipbot
340 write (iout,*) "SHIELD MODE",shield_mode
341 if (TUBElog.gt.0) then
342 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
343 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
344 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
345 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
346 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
347 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
348 buftubebot=bordtubebot+tubebufthick
349 buftubetop=bordtubetop-tubebufthick
351 if (me.eq.king .or. .not.out1file )
352 & write (iout,*) "DISTCHAINMAX",distchainmax
354 if(me.eq.king.or..not.out1file)
355 & write (iout,'(2a)') diagmeth(kdiag),
356 & ' routine used to diagonalize matrices.'
359 c--------------------------------------------------------------------------
360 subroutine read_REMDpar
366 include 'COMMON.IOUNITS'
367 include 'COMMON.TIME1'
370 include 'COMMON.LANGEVIN'
373 include 'COMMON.LANGEVIN.lang0.5diag'
375 include 'COMMON.LANGEVIN.lang0'
378 include 'COMMON.INTERACT'
379 include 'COMMON.NAMES'
381 include 'COMMON.REMD'
382 include 'COMMON.CONTROL'
383 include 'COMMON.SETUP'
385 character*320 controlcard
386 character*3200 controlcard1
387 integer iremd_m_total,i
389 if(me.eq.king.or..not.out1file)
390 & write (iout,*) "REMD setup"
392 call card_concat(controlcard)
393 call readi(controlcard,"NREP",nrep,3)
394 call readi(controlcard,"NSTEX",nstex,1000)
395 call reada(controlcard,"RETMIN",retmin,10.0d0)
396 call reada(controlcard,"RETMAX",retmax,1000.0d0)
397 mremdsync=(index(controlcard,'SYNC').gt.0)
398 call readi(controlcard,"NSYN",i_sync_step,100)
399 restart1file=(index(controlcard,'REST1FILE').gt.0)
400 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
401 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
402 if(max_cache_traj_use.gt.max_cache_traj)
403 & max_cache_traj_use=max_cache_traj
404 if(me.eq.king.or..not.out1file) then
405 cd if (traj1file) then
406 crc caching is in testing - NTWX is not ignored
407 cd write (iout,*) "NTWX value is ignored"
408 cd write (iout,*) " trajectory is stored to one file by master"
409 cd write (iout,*) " before exchange at NSTEX intervals"
411 write (iout,*) "NREP= ",nrep
412 write (iout,*) "NSTEX= ",nstex
413 write (iout,*) "SYNC= ",mremdsync
414 write (iout,*) "NSYN= ",i_sync_step
415 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
418 if (index(controlcard,'TLIST').gt.0) then
420 call card_concat(controlcard1)
421 read(controlcard1,*) (remd_t(i),i=1,nrep)
422 if(me.eq.king.or..not.out1file)
423 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
426 if (index(controlcard,'MLIST').gt.0) then
428 call card_concat(controlcard1)
429 read(controlcard1,*) (remd_m(i),i=1,nrep)
430 if(me.eq.king.or..not.out1file) then
431 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
434 iremd_m_total=iremd_m_total+remd_m(i)
436 write (iout,*) 'Total number of replicas ',iremd_m_total
441 & "Adaptive (PMF-biased) umbrella sampling will be run"
444 if(me.eq.king.or..not.out1file)
445 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
448 c--------------------------------------------------------------------------
449 subroutine read_MDpar
455 include 'COMMON.IOUNITS'
456 include 'COMMON.TIME1'
458 include 'COMMON.QRESTR'
460 include 'COMMON.LANGEVIN'
463 include 'COMMON.LANGEVIN.lang0.5diag'
465 include 'COMMON.LANGEVIN.lang0'
468 include 'COMMON.INTERACT'
469 include 'COMMON.NAMES'
471 include 'COMMON.SETUP'
472 include 'COMMON.CONTROL'
473 include 'COMMON.SPLITELE'
474 include 'COMMON.FFIELD'
476 character*320 controlcard
480 call card_concat(controlcard)
481 call readi(controlcard,"NSTEP",n_timestep,1000000)
482 call readi(controlcard,"NTWE",ntwe,100)
483 call readi(controlcard,"NTWX",ntwx,1000)
484 call reada(controlcard,"DT",d_time,1.0d-1)
485 call reada(controlcard,"DVMAX",dvmax,2.0d1)
486 call reada(controlcard,"DAMAX",damax,1.0d1)
487 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
488 call readi(controlcard,"LANG",lang,0)
489 RESPA = index(controlcard,"RESPA") .gt. 0
490 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
491 ntime_split0=ntime_split
492 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
493 ntime_split0=ntime_split
494 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
495 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
496 rest = index(controlcard,"REST").gt.0
497 tbf = index(controlcard,"TBF").gt.0
498 usampl = index(controlcard,"USAMPL").gt.0
499 scale_umb = index(controlcard,"SCALE_UMB").gt.0
500 adaptive = index(controlcard,"ADAPTIVE").gt.0
501 mdpdb = index(controlcard,"MDPDB").gt.0
502 call reada(controlcard,"T_BATH",t_bath,300.0d0)
503 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
504 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
505 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
506 if (count_reset_moment.eq.0) count_reset_moment=1000000000
507 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
508 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
509 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
510 if (count_reset_vel.eq.0) count_reset_vel=1000000000
511 large = index(controlcard,"LARGE").gt.0
512 print_compon = index(controlcard,"PRINT_COMPON").gt.0
513 rattle = index(controlcard,"RATTLE").gt.0
514 preminim = index(controlcard,"PREMINIM").gt.0
515 if (iranconf.gt.0 .or. indpdb.gt.0 .or. start_from_model) then
519 write (iout,*) "PREMINIM ",preminim
521 if (index(controlcard,'CART').gt.0) dccart=.true.
522 write (iout,*) "dccart ",dccart
523 write (iout,*) "read_minim ",index(controlcard,'READ_MINIM_PAR')
524 if (index(controlcard,'READ_MINIM_PAR').gt.0) call read_minim
526 c if performing umbrella sampling, fragments constrained are read from the fragment file
529 write (iout,*) "Umbrella sampling will be run"
530 if (scale_umb.and.adaptive) then
531 write (iout,*) "ADAPTIVE and SCALE_UMB are mutually exclusive"
532 write (iout,*) "Select one of those and re-run the job."
535 if (scale_umb) write (iout,*)
536 &"Umbrella-restraint force constants will be scaled by temperature"
540 if(me.eq.king.or..not.out1file) then
542 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
544 write (iout,'(a)') "The units are:"
545 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
546 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
547 & " acceleration: angstrom/(48.9 fs)**2"
548 write (iout,'(a)') "energy: kcal/mol, temperature: K"
550 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
551 write (iout,'(a60,f10.5,a)')
552 & "Initial time step of numerical integration:",d_time,
554 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
555 write (iout,'(a60,f10.5,a)') "Cutoff on interactions",r_cut_int,
557 write(iout,'(a60,i5)')"Frequency of updating interaction list",
560 write (iout,'(2a,i4,a)')
561 & "A-MTS algorithm used; initial time step for fast-varying",
562 & " short-range forces split into",ntime_split," steps."
563 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
564 & r_cut_respa," lambda",rlamb
566 write (iout,'(2a,f10.5)')
567 & "Maximum acceleration threshold to reduce the time step",
568 & "/increase split number:",damax
569 write (iout,'(2a,f10.5)')
570 & "Maximum predicted energy drift to reduce the timestep",
571 & "/increase split number:",edriftmax
572 write (iout,'(a60,f10.5)')
573 & "Maximum velocity threshold to reduce velocities:",dvmax
574 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
575 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
576 if (rattle) write (iout,'(a60)')
577 & "Rattle algorithm used to constrain the virtual bonds"
581 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
582 call reada(controlcard,"RWAT",rwat,1.4d0)
583 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
584 surfarea=index(controlcard,"SURFAREA").gt.0
585 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
586 if(me.eq.king.or..not.out1file)then
587 write (iout,'(/a,$)') "Langevin dynamics calculation"
590 & " with direct integration of Langevin equations"
591 else if (lang.eq.2) then
592 write (iout,'(a/)') " with TINKER stochasic MD integrator"
593 else if (lang.eq.3) then
594 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
595 else if (lang.eq.4) then
596 write (iout,'(a/)') " in overdamped mode"
598 write (iout,'(//a,i5)')
599 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
602 write (iout,'(a60,f10.5)') "Temperature:",t_bath
603 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
604 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
605 write (iout,'(a60,f10.5)')
606 & "Scaling factor of the friction forces:",scal_fric
607 if (surfarea) write (iout,'(2a,i10,a)')
608 & "Friction coefficients will be scaled by solvent-accessible",
609 & " surface area every",reset_fricmat," steps."
611 c Calculate friction coefficients and bounds of stochastic forces
612 eta=6*pi*cPoise*etawat
613 if(me.eq.king.or..not.out1file)
614 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
616 gamp=scal_fric*(pstok+rwat)*eta
617 stdfp=dsqrt(2*Rb*t_bath/d_time)
619 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
620 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
622 if(me.eq.king.or..not.out1file)then
623 write (iout,'(/2a/)')
624 & "Radii of site types and friction coefficients and std's of",
625 & " stochastic forces of fully exposed sites"
626 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
628 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
629 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
633 if(me.eq.king.or..not.out1file)then
634 write (iout,'(a)') "Berendsen bath calculation"
635 write (iout,'(a60,f10.5)') "Temperature:",t_bath
636 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
638 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
639 & count_reset_moment," steps"
641 & write (iout,'(a,i10,a)')
642 & "Velocities will be reset at random every",count_reset_vel,
646 if(me.eq.king.or..not.out1file)
647 & write (iout,'(a31)') "Microcanonical mode calculation"
649 if(me.eq.king.or..not.out1file)then
650 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
652 write(iout,*) "MD running with constraints."
653 write(iout,*) "Equilibration time ", eq_time, " mtus."
654 write(iout,*) "Constraining ", nfrag," fragments."
655 write(iout,*) "Length of each fragment, weight and q0:"
657 write (iout,*) "Set of restraints #",iset
659 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
660 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
662 write(iout,*) "constraints between ", npair, "fragments."
663 write(iout,*) "constraint pairs, weights and q0:"
665 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
666 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
668 write(iout,*) "angle constraints within ", nfrag_back,
669 & "backbone fragments."
671 write(iout,*) "fragment, weights, q0:"
673 write(iout,'(2i5,3(f8.1,f8.2))') ifrag_back(1,i,iset),
674 & ifrag_back(2,i,iset),
675 & wfrag_back(1,i,iset),qin_back(1,i,iset),
676 & wfrag_back(2,i,iset),qin_back(2,i,iset),
677 & wfrag_back(3,i,iset),qin_back(3,i,iset)
680 write(iout,*) "fragment, weights:"
682 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
683 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
684 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
688 iset=mod(kolor,nset)+1
691 if(me.eq.king.or..not.out1file)
692 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
695 c------------------------------------------------------------------------------
698 C Read molecular data.
704 integer error_msg,ierror,ierr,ierrcode
706 include 'COMMON.IOUNITS'
709 include 'COMMON.INTERACT'
710 include 'COMMON.LOCAL'
711 include 'COMMON.NAMES'
712 include 'COMMON.CHAIN'
713 include 'COMMON.FFIELD'
714 include 'COMMON.SBRIDGE'
715 include 'COMMON.HEADER'
716 include 'COMMON.CONTROL'
717 include 'COMMON.SAXS'
718 include 'COMMON.DBASE'
719 include 'COMMON.THREAD'
720 include 'COMMON.CONTACTS'
721 include 'COMMON.TORCNSTR'
722 include 'COMMON.TIME1'
723 include 'COMMON.BOUNDS'
725 include 'COMMON.SETUP'
726 include 'COMMON.SHIELD'
727 character*4 sequence(maxres)
729 double precision x(maxvar)
730 character*256 pdbfile
731 character*400 weightcard
732 character*80 weightcard_t,ucase
733 integer itype_pdb(maxres)
734 common /pizda/ itype_pdb
735 logical seq_comp,fail
736 double precision energia(0:n_ene)
737 double precision secprob(3,maxdih_constr)
739 double precision phihel,phibet,sigmahel,sigmabet
740 integer iti,nsi,maxsi
744 integer i,j,ii,k,l,itrial,itmp,i1,i2,it1,it2,nres_temp
745 double precision sumv
747 C Read PDB structure if applicable
749 if (indpdb.gt.0 .or. pdbref) then
750 read(inp,'(a)') pdbfile
751 if(me.eq.king.or..not.out1file)
752 & write (iout,'(2a)') 'PDB data will be read from file ',
753 & pdbfile(:ilen(pdbfile))
754 open(ipdbin,file=pdbfile,status='old',err=33)
756 33 write (iout,'(a)') 'Error opening PDB file.'
767 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
768 & (crefjlee(j,i+nres),j=1,3)
771 if(me.eq.king.or..not.out1file)
772 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
773 & ' nstart_sup=',nstart_sup
775 itype_pdb(i)=itype(i)
779 nct=nstart_sup+nsup-1
780 call contact(.false.,ncont_ref,icont_ref,co)
783 if(me.eq.king.or..not.out1file)
784 & write(iout,*)'Adding sidechains'
788 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
791 do while (fail.and.nsi.le.maxsi)
792 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
795 if(fail) write(iout,*)'Adding sidechain failed for res ',
796 & i,' after ',nsi,' trials'
801 if (indpdb.eq.0) then
802 C Read sequence if not taken from the pdb file.
804 c print *,'nres=',nres
805 if (iscode.gt.0) then
806 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
808 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
810 C Convert sequence to numeric code
812 itype(i)=rescode(i,sequence(i),iscode)
816 c print '(20i4)',(itype(i),i=1,nres)
819 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
821 if (itype(i).eq.ntyp1) then
825 else if (iabs(itype(i+1)).ne.20) then
827 else if (iabs(itype(i)).ne.20) then
834 if(me.eq.king.or..not.out1file)then
835 write (iout,*) "ITEL"
837 write (iout,*) i,itype(i),itel(i)
839 c print *,'Call Read_Bridge.'
843 cd print *,'NNT=',NNT,' NCT=',NCT
844 call seq2chains(nres,itype,nchain,chain_length,chain_border,
847 chain_border1(2,1)=chain_border(2,1)+1
849 chain_border1(1,i)=chain_border(1,i)-1
850 chain_border1(2,i)=chain_border(2,i)+1
852 if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
853 chain_border1(2,nchain)=nres
854 write(iout,*) "nres",nres," nchain",nchain
856 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),
857 & chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
859 call chain_symmetry(nchain,nres,itype,chain_border,
860 & chain_length,npermchain,tabpermchain)
862 c write(iout,*) i,(tperm(ireschain(i),ii,tabpermchain),
865 write(iout,*) "residue permutations"
867 write(iout,*) i,(iperm(i,ii),ii=1,npermchain)
870 if (itype(1).eq.ntyp1) nnt=2
871 if (itype(nres).eq.ntyp1) nct=nct-1
872 write (iout,*) "nnt",nnt," nct",nct
875 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
876 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
878 c print*, 'init_dfa_vars finished!'
880 c print*, 'read_dfa_info finished!'
884 if(me.eq.king.or..not.out1file)
885 & write (iout,'(a,i3)') 'nsup=',nsup
887 if (nsup.le.(nct-nnt+1)) then
888 do i=0,nct-nnt+1-nsup
889 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
895 & 'Error - sequences to be superposed do not match.'
898 do i=0,nsup-(nct-nnt+1)
899 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
901 nstart_sup=nstart_sup+i
907 & 'Error - sequences to be superposed do not match.'
910 if (nsup.eq.0) nsup=nct-nnt
911 if (nstart_sup.eq.0) nstart_sup=nnt
912 if (nstart_seq.eq.0) nstart_seq=nnt
913 if(me.eq.king.or..not.out1file)
914 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
915 & ' nstart_seq=',nstart_seq
918 C 8/13/98 Set limits to generating the dihedral angles
923 read (inp,*) ndih_constr
924 write (iout,*) "ndih_constr",ndih_constr
925 if (ndih_constr.gt.0) then
928 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
931 if(me.eq.king.or..not.out1file)then
934 & 'There are',ndih_constr,' restraints on gamma angles.'
936 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
943 phi0(i)=deg2rad*phi0(i)
944 drange(i)=deg2rad*drange(i)
946 C if(me.eq.king.or..not.out1file)
947 C & write (iout,*) 'FTORS',ftors
950 phibound(1,ii) = phi0(i)-drange(i)
951 phibound(2,ii) = phi0(i)+drange(i)
954 if (me.eq.king .or. .not.out1file) then
956 write (iout,'(a)') 'Boundaries in gamma angle sampling:'
958 write (iout,'(a3,i5,2f10.1)')
959 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
964 else if (ndih_constr.lt.0) then
966 call card_concat(weightcard)
967 call reada(weightcard,"PHIHEL",phihel,50.0D0)
968 call reada(weightcard,"PHIBET",phibet,180.0D0)
969 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
970 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
971 call reada(weightcard,"WDIHC",wdihc,0.591D0)
972 write (iout,*) "Weight of dihedral angle restraints",wdihc
973 read(inp,'(9x,3f7.3)')
974 & (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
975 write (iout,*) "The secprob array"
977 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
981 if (itype(i-3).ne.ntyp1 .and. itype(i-2).ne.ntyp1
982 & .and. itype(i-1).ne.ntyp1 .and. itype(i).ne.ntyp1) then
983 ndih_constr=ndih_constr+1
984 idih_constr(ndih_constr)=i
987 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
988 sumv=sumv+vpsipred(j,ndih_constr)
991 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
993 phibound(1,ndih_constr)=phihel*deg2rad
994 phibound(2,ndih_constr)=phibet*deg2rad
995 sdihed(1,ndih_constr)=sigmahel*deg2rad
996 sdihed(2,ndih_constr)=sigmabet*deg2rad
1000 if(me.eq.king.or..not.out1file)then
1003 & 'There are',ndih_constr,
1004 & ' bimodal restraints on gamma angles.'
1006 write(iout,'(i5,1x,a4,i5,1h-,a4,i5,4f8.3,3f10.5)') i,
1007 & restyp(itype(idih_constr(i)-2)),idih_constr(i)-2,
1008 & restyp(itype(idih_constr(i)-1)),idih_constr(i)-1,
1009 & phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,
1010 & phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,
1011 & (vpsipred(j,i),j=1,3)
1018 C first setting the theta boundaries to 0 to pi
1019 C this mean that there is no energy penalty for any angle occuring this can be applied
1020 C for generate random conformation but is not implemented in this way
1023 C thetabound(2,i)=pi
1025 C begin reading theta constrains this is quartic constrains allowing to
1026 C have smooth second derivative
1027 if (with_theta_constr) then
1028 C with_theta_constr is keyword allowing for occurance of theta constrains
1029 read (inp,*) ntheta_constr
1030 C ntheta_constr is the number of theta constrains
1031 if (ntheta_constr.gt.0) then
1032 C read (inp,*) ftors
1033 read (inp,*) (itheta_constr(i),theta_constr0(i),
1034 & theta_drange(i),for_thet_constr(i),
1035 & i=1,ntheta_constr)
1036 C the above code reads from 1 to ntheta_constr
1037 C itheta_constr(i) residue i for which is theta_constr
1038 C theta_constr0 the global minimum value
1039 C theta_drange is range for which there is no energy penalty
1040 C for_thet_constr is the force constant for quartic energy penalty
1042 if(me.eq.king.or..not.out1file)then
1044 & 'There are',ntheta_constr,' constraints on phi angles.'
1045 do i=1,ntheta_constr
1046 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
1048 & for_thet_constr(i)
1051 do i=1,ntheta_constr
1052 theta_constr0(i)=deg2rad*theta_constr0(i)
1053 theta_drange(i)=deg2rad*theta_drange(i)
1055 C if(me.eq.king.or..not.out1file)
1056 C & write (iout,*) 'FTORS',ftors
1057 C do i=1,ntheta_constr
1058 C ii = itheta_constr(i)
1059 C thetabound(1,ii) = phi0(i)-drange(i)
1060 C thetabound(2,ii) = phi0(i)+drange(i)
1062 endif ! ntheta_constr.gt.0
1063 endif! with_theta_constr
1065 C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
1066 C write (iout,*) "with_dihed_constr ",with_dihed_constr
1067 c--- Zscore rms -------
1068 if (nz_start.eq.0) nz_start=nnt
1069 if (nz_end.eq.0 .and. nsup.gt.0) then
1071 else if (nz_end.eq.0) then
1074 if(me.eq.king.or..not.out1file)then
1075 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1076 write (iout,*) 'IZ_SC=',iz_sc
1078 c----------------------
1081 if (.not.pdbref) then
1082 call read_angles(inp,*38)
1085 38 write (iout,'(a)') 'Error reading reference structure.'
1087 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1088 stop 'Error reading reference structure'
1090 39 call chainbuild_extconf
1092 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1101 call contact(.true.,ncont_ref,icont_ref,co)
1105 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1107 if (constr_dist.gt.0) call read_dist_constr
1108 c write (iout,*) "After read_dist_constr nhpb",nhpb
1109 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1111 call NMRpeak_partition
1112 if(me.eq.king.or..not.out1file)write (iout,*) 'Contact order:',co
1114 if(me.eq.king.or..not.out1file)
1115 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1118 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1120 if(me.eq.king.or..not.out1file)
1121 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1122 & icont_ref(1,i),' ',
1123 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1126 write (iout,*) "calling read_saxs_consrtr",nsaxs
1127 if (nsaxs.gt.0) call read_saxs_constr
1129 if (constr_homology.gt.0) then
1130 call read_constr_homology
1131 if (indpdb.gt.0 .or. pdbref) then
1134 c(j,i)=crefjlee(j,i)
1135 cref(j,i)=crefjlee(j,i)
1140 write (iout,*) "sc_loc_geom: Array C"
1142 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1143 & (c(j,i+nres),j=1,3)
1145 write (iout,*) "Array Cref"
1147 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1148 & (cref(j,i+nres),j=1,3)
1151 call int_from_cart1(.false.)
1152 call sc_loc_geom(.false.)
1154 thetaref(i)=theta(i)
1159 dc(j,i)=c(j,i+1)-c(j,i)
1160 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1165 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1166 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1171 if (start_from_model) then
1174 read(inp,'(a)',end=332,err=332) pdbfile
1175 if (me.eq.king .or. .not. out1file)
1176 & write (iout,'(a,5x,a)') 'Opening PDB file',
1177 & pdbfile(:ilen(pdbfile))
1178 open(ipdbin,file=pdbfile,status='old',err=336)
1180 336 write (iout,'(a,5x,a)') 'Error opening PDB file',
1181 & pdbfile(:ilen(pdbfile))
1189 if (nres.ge.nres_temp) then
1190 nmodel_start=nmodel_start+1
1191 pdbfiles_chomo(nmodel_start)=pdbfile
1194 chomo(j,i,nmodel_start)=c(j,i)
1198 if (me.eq.king .or. .not. out1file)
1199 & write (iout,'(a,2i5,1x,a)')
1200 & "Different number of residues",nres_temp,nres,
1206 if (nmodel_start.eq.0) then
1207 if (me.eq.king .or. .not. out1file)
1208 & write (iout,'(a)')
1209 & "No valid starting model found START_FROM_MODELS is OFF"
1210 start_from_model=.false.
1212 write (iout,*) "nmodel_start",nmodel_start
1218 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1219 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1220 & modecalc.ne.10) then
1221 C If input structure hasn't been supplied from the PDB file read or generate
1223 if (iranconf.eq.0 .and. .not. extconf .and. .not.
1224 & start_from_model) then
1225 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1226 & write (iout,'(a)') 'Initial geometry will be read in.'
1228 read(inp,'(8f10.5)',end=36,err=36)
1229 & ((c(l,k),l=1,3),k=1,nres),
1230 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1231 c write (iout,*) "Exit READ_CART"
1232 c write (iout,'(8f10.5)')
1233 c & ((c(l,k),l=1,3),k=1,nres),
1234 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1238 dc(j,i)=c(j,i+1)-c(j,i)
1239 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1243 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1245 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1246 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1251 c dc_norm(j,i+nres)=0.0d0
1255 call int_from_cart1(.true.)
1256 write (iout,*) "Finish INT_TO_CART"
1257 c write (iout,*) "DC"
1259 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1260 c & (dc(j,i+nres),j=1,3)
1266 call read_angles(inp,*36)
1268 call chainbuild_extconf
1271 36 write (iout,'(a)') 'Error reading angle file.'
1273 call mpi_finalize( MPI_COMM_WORLD,IERR )
1275 stop 'Error reading angle file.'
1277 else if (extconf) then
1278 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1279 & write (iout,'(a)') 'Extended chain initial geometry.'
1281 theta(i)=90d0*deg2rad
1284 phi(i)=180d0*deg2rad
1287 alph(i)=110d0*deg2rad
1290 omeg(i)=-120d0*deg2rad
1291 if (itype(i).le.0) omeg(i)=-omeg(i)
1294 call chainbuild_extconf
1296 if(me.eq.king.or..not.out1file)
1297 & write (iout,'(a)') 'Random-generated initial geometry.'
1300 if (me.eq.king .or. fg_rank.eq.0 .and. (
1301 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1305 call gen_rand_conf(itmp,*30)
1307 30 write (iout,*) 'Failed to generate random conformation',
1308 & ', itrial=',itrial
1309 write (*,*) 'Processor:',me,
1310 & ' Failed to generate random conformation',
1320 write (iout,'(a,i3,a)') 'Processor:',me,
1321 & ' error in generating random conformation.'
1322 write (*,'(a,i3,a)') 'Processor:',me,
1323 & ' error in generating random conformation.'
1326 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1331 & ' error in generating random conformation.'
1336 elseif (modecalc.eq.4) then
1337 read (inp,'(a)') intinname
1338 open (intin,file=intinname,status='old',err=333)
1339 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1340 & write (iout,'(a)') 'intinname',intinname
1341 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1343 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1345 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1347 stop 'Error opening angle file.'
1351 C Generate distance constraints, if the PDB structure is to be regularized.
1352 if (nthread.gt.0) then
1353 call read_threadbase
1356 if (me.eq.king .or. .not. out1file)
1358 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1359 write (iout,'(/a,i3,a)')
1360 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1361 write (iout,'(20i4)') (iss(i),i=1,ns)
1363 write(iout,*)"Running with dynamic disulfide-bond formation"
1365 write (iout,'(/a/)') 'Pre-formed links are:'
1371 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1372 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1378 if (ns.gt.0.and.dyn_ss) then
1382 forcon(i-nss)=forcon(i)
1389 dyn_ss_mask(iss(i))=.true.
1394 c write (iout,*) "DC"
1396 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1397 c & (dc(j,i+nres),j=1,3)
1399 if (i2ndstr.gt.0) call secstrp2dihc
1400 c call geom_to_var(nvar,x)
1401 c call etotal(energia(0))
1402 c call enerprint(energia(0))
1403 c call briefout(0,etot)
1405 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1406 cd write (iout,'(a)') 'Variable list:'
1407 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1409 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1410 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1411 & 'Processor',myrank,': end reading molecular data.'
1416 c--------------------------------------------------------------------------
1417 logical function seq_comp(itypea,itypeb,length)
1419 integer length,itypea(length),itypeb(length)
1422 if (itypea(i).ne.itypeb(i)) then
1430 c-----------------------------------------------------------------------------
1431 subroutine read_bridge
1432 C Read information about disulfide bridges.
1434 include 'DIMENSIONS'
1439 include 'COMMON.IOUNITS'
1440 include 'COMMON.GEO'
1441 include 'COMMON.VAR'
1442 include 'COMMON.INTERACT'
1443 include 'COMMON.LOCAL'
1444 include 'COMMON.NAMES'
1445 include 'COMMON.CHAIN'
1446 include 'COMMON.FFIELD'
1447 include 'COMMON.SBRIDGE'
1448 include 'COMMON.HEADER'
1449 include 'COMMON.CONTROL'
1450 include 'COMMON.DBASE'
1451 include 'COMMON.THREAD'
1452 include 'COMMON.TIME1'
1453 include 'COMMON.SETUP'
1455 C Read bridging residues.
1456 read (inp,*) ns,(iss(i),i=1,ns)
1457 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
1464 if(me.eq.king.or..not.out1file)
1465 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1466 C Check whether the specified bridging residues are cystines.
1468 if (itype(iss(i)).ne.1) then
1469 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1470 & 'Do you REALLY think that the residue ',
1471 & restyp(itype(iss(i))),i,
1472 & ' can form a disulfide bridge?!!!'
1473 write (*,'(2a,i3,a)')
1474 & 'Do you REALLY think that the residue ',
1475 & restyp(itype(iss(i))),i,
1476 & ' can form a disulfide bridge?!!!'
1478 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1483 C Read preformed bridges.
1485 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1487 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1490 C Check if the residues involved in bridges are in the specified list of
1491 C bridging residues.
1494 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1495 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1496 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1497 & ' contains residues present in other pairs.'
1498 write (*,'(a,i3,a)') 'Disulfide pair',i,
1499 & ' contains residues present in other pairs.'
1501 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1507 if (ihpb(i).eq.iss(j)) goto 10
1509 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1512 if (jhpb(i).eq.iss(j)) goto 20
1514 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1520 ihpb(i)=ihpb(i)+nres
1521 jhpb(i)=jhpb(i)+nres
1527 c----------------------------------------------------------------------------
1528 subroutine read_x(kanal,*)
1530 include 'DIMENSIONS'
1531 include 'COMMON.GEO'
1532 include 'COMMON.VAR'
1533 include 'COMMON.CHAIN'
1534 include 'COMMON.IOUNITS'
1535 include 'COMMON.CONTROL'
1536 include 'COMMON.LOCAL'
1537 include 'COMMON.INTERACT'
1538 integer i,j,k,l,kanal
1539 c Read coordinates from input
1541 read(kanal,'(8f10.5)',end=10,err=10)
1542 & ((c(l,k),l=1,3),k=1,nres),
1543 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1546 c(j,2*nres)=c(j,nres)
1548 call int_from_cart1(.false.)
1551 dc(j,i)=c(j,i+1)-c(j,i)
1552 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1556 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1558 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1559 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1567 c----------------------------------------------------------------------------
1568 subroutine read_threadbase
1570 include 'DIMENSIONS'
1571 include 'COMMON.IOUNITS'
1572 include 'COMMON.GEO'
1573 include 'COMMON.VAR'
1574 include 'COMMON.INTERACT'
1575 include 'COMMON.LOCAL'
1576 include 'COMMON.NAMES'
1577 include 'COMMON.CHAIN'
1578 include 'COMMON.FFIELD'
1579 include 'COMMON.SBRIDGE'
1580 include 'COMMON.HEADER'
1581 include 'COMMON.CONTROL'
1582 include 'COMMON.DBASE'
1583 include 'COMMON.THREAD'
1584 include 'COMMON.TIME1'
1586 double precision dist
1587 C Read pattern database for threading.
1588 read (icbase,*) nseq
1590 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1591 & nres_base(2,i),nres_base(3,i)
1592 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1594 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1595 c & nres_base(2,i),nres_base(3,i)
1596 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1600 if (weidis.eq.0.0D0) weidis=0.1D0
1609 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1610 write (iout,'(a,i5)') 'nexcl: ',nexcl
1611 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1614 c------------------------------------------------------------------------------
1615 subroutine setup_var
1617 include 'DIMENSIONS'
1618 include 'COMMON.IOUNITS'
1619 include 'COMMON.GEO'
1620 include 'COMMON.VAR'
1621 include 'COMMON.INTERACT'
1622 include 'COMMON.LOCAL'
1623 include 'COMMON.NAMES'
1624 include 'COMMON.CHAIN'
1625 include 'COMMON.FFIELD'
1626 include 'COMMON.SBRIDGE'
1627 include 'COMMON.HEADER'
1628 include 'COMMON.CONTROL'
1629 include 'COMMON.DBASE'
1630 include 'COMMON.THREAD'
1631 include 'COMMON.TIME1'
1633 C Set up variable list.
1639 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1641 ialph(i,1)=nvar+nside
1645 if (indphi.gt.0) then
1647 else if (indback.gt.0) then
1654 c----------------------------------------------------------------------------
1655 subroutine gen_dist_constr
1656 C Generate CA distance constraints.
1658 include 'DIMENSIONS'
1659 include 'COMMON.IOUNITS'
1660 include 'COMMON.GEO'
1661 include 'COMMON.VAR'
1662 include 'COMMON.INTERACT'
1663 include 'COMMON.LOCAL'
1664 include 'COMMON.NAMES'
1665 include 'COMMON.CHAIN'
1666 include 'COMMON.FFIELD'
1667 include 'COMMON.SBRIDGE'
1668 include 'COMMON.HEADER'
1669 include 'COMMON.CONTROL'
1670 include 'COMMON.DBASE'
1671 include 'COMMON.THREAD'
1672 include 'COMMON.SPLITELE'
1673 include 'COMMON.TIME1'
1674 integer i,j,itype_pdb(maxres)
1675 common /pizda/ itype_pdb
1677 double precision dist
1679 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1680 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1681 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1683 do i=nstart_sup,nstart_sup+nsup-1
1684 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1685 cd & ' seq_pdb', restyp(itype_pdb(i))
1686 do j=i+2,nstart_sup+nsup-1
1687 c 5/24/2020 Adam: Cutoff included to reduce array size
1689 if (dd.gt.r_cut_int) cycle
1691 ihpb(nhpb)=i+nstart_seq-nstart_sup
1692 jhpb(nhpb)=j+nstart_seq-nstart_sup
1697 cd write (iout,'(a)') 'Distance constraints:'
1702 cd if (ii.gt.nres) then
1707 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1708 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1709 cd & dhpb(i),forcon(i)
1713 c----------------------------------------------------------------------------
1716 include 'DIMENSIONS'
1717 include 'COMMON.MAP'
1718 include 'COMMON.IOUNITS'
1720 character*3 angid(4) /'THE','PHI','ALP','OME'/
1721 character*80 mapcard,ucase
1723 read (inp,'(a)') mapcard
1724 mapcard=ucase(mapcard)
1725 if (index(mapcard,'PHI').gt.0) then
1727 else if (index(mapcard,'THE').gt.0) then
1729 else if (index(mapcard,'ALP').gt.0) then
1731 else if (index(mapcard,'OME').gt.0) then
1734 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1735 stop 'Error - illegal variable spec in MAP card.'
1737 call readi (mapcard,'RES1',res1(imap),0)
1738 call readi (mapcard,'RES2',res2(imap),0)
1739 if (res1(imap).eq.0) then
1740 res1(imap)=res2(imap)
1741 else if (res2(imap).eq.0) then
1742 res2(imap)=res1(imap)
1744 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1746 & 'Error - illegal definition of variable group in MAP.'
1747 stop 'Error - illegal definition of variable group in MAP.'
1749 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1750 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1751 call readi(mapcard,'NSTEP',nstep(imap),0)
1752 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1754 & 'Illegal boundary and/or step size specification in MAP.'
1755 stop 'Illegal boundary and/or step size specification in MAP.'
1760 c----------------------------------------------------------------------------
1763 include 'DIMENSIONS'
1764 include 'COMMON.IOUNITS'
1765 include 'COMMON.GEO'
1766 include 'COMMON.CSA'
1767 include 'COMMON.BANK'
1768 include 'COMMON.CONTROL'
1770 character*620 mcmcard
1771 call card_concat(mcmcard)
1773 call readi(mcmcard,'NCONF',nconf,50)
1774 call readi(mcmcard,'NADD',nadd,0)
1775 call readi(mcmcard,'JSTART',jstart,1)
1776 call readi(mcmcard,'JEND',jend,1)
1777 call readi(mcmcard,'NSTMAX',nstmax,500000)
1778 call readi(mcmcard,'N0',n0,1)
1779 call readi(mcmcard,'N1',n1,6)
1780 call readi(mcmcard,'N2',n2,4)
1781 call readi(mcmcard,'N3',n3,0)
1782 call readi(mcmcard,'N4',n4,0)
1783 call readi(mcmcard,'N5',n5,0)
1784 call readi(mcmcard,'N6',n6,10)
1785 call readi(mcmcard,'N7',n7,0)
1786 call readi(mcmcard,'N8',n8,0)
1787 call readi(mcmcard,'N9',n9,0)
1788 call readi(mcmcard,'N14',n14,0)
1789 call readi(mcmcard,'N15',n15,0)
1790 call readi(mcmcard,'N16',n16,0)
1791 call readi(mcmcard,'N17',n17,0)
1792 call readi(mcmcard,'N18',n18,0)
1794 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1796 call readi(mcmcard,'NDIFF',ndiff,2)
1797 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1798 call readi(mcmcard,'IS1',is1,1)
1799 call readi(mcmcard,'IS2',is2,8)
1800 call readi(mcmcard,'NRAN0',nran0,4)
1801 call readi(mcmcard,'NRAN1',nran1,2)
1802 call readi(mcmcard,'IRR',irr,1)
1803 call readi(mcmcard,'NSEED',nseed,20)
1804 call readi(mcmcard,'NTOTAL',ntotal,10000)
1805 call reada(mcmcard,'CUT1',cut1,2.0d0)
1806 call reada(mcmcard,'CUT2',cut2,5.0d0)
1807 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1808 call readi(mcmcard,'ICMAX',icmax,3)
1809 call readi(mcmcard,'IRESTART',irestart,0)
1810 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1813 call reada(mcmcard,'DELE',dele,20.0d0)
1814 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1815 call readi(mcmcard,'IREF',iref,0)
1816 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1817 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1818 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1819 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1820 write (iout,*) "NCONF_IN",nconf_in
1823 c----------------------------------------------------------------------------
1826 include 'DIMENSIONS'
1827 include 'COMMON.MCM'
1828 include 'COMMON.MCE'
1829 include 'COMMON.IOUNITS'
1831 character*320 mcmcard
1833 call card_concat(mcmcard)
1834 call readi(mcmcard,'MAXACC',maxacc,100)
1835 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1836 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1837 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1838 call readi(mcmcard,'MAXREPM',maxrepm,200)
1839 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1840 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1841 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1842 call reada(mcmcard,'E_UP',e_up,5.0D0)
1843 call reada(mcmcard,'DELTE',delte,0.1D0)
1844 call readi(mcmcard,'NSWEEP',nsweep,5)
1845 call readi(mcmcard,'NSTEPH',nsteph,0)
1846 call readi(mcmcard,'NSTEPC',nstepc,0)
1847 call reada(mcmcard,'TMIN',tmin,298.0D0)
1848 call reada(mcmcard,'TMAX',tmax,298.0D0)
1849 call readi(mcmcard,'NWINDOW',nwindow,0)
1850 call readi(mcmcard,'PRINT_MC',print_mc,0)
1851 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1852 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1853 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1854 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1855 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1856 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1857 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1858 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1859 if (nwindow.gt.0) then
1860 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1862 winlen(i)=winend(i)-winstart(i)+1
1865 if (tmax.lt.tmin) tmax=tmin
1866 if (tmax.eq.tmin) then
1870 if (nstepc.gt.0 .and. nsteph.gt.0) then
1871 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1872 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1874 C Probabilities of different move types
1875 sumpro_type(0)=0.0D0
1876 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1877 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1878 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1879 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1880 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1881 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1882 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1884 print *,'i',i,' sumprotype',sumpro_type(i)
1885 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1886 print *,'i',i,' sumprotype',sumpro_type(i)
1890 c----------------------------------------------------------------------------
1891 subroutine read_minim
1893 include 'DIMENSIONS'
1894 include 'COMMON.MINIM'
1895 include 'COMMON.IOUNITS'
1897 character*320 minimcard
1898 call card_concat(minimcard)
1899 call readi(minimcard,'MAXMIN',maxmin,2000)
1900 call readi(minimcard,'MAXFUN',maxfun,5000)
1901 call readi(minimcard,'MINMIN',minmin,maxmin)
1902 call readi(minimcard,'MINFUN',minfun,maxmin)
1903 call reada(minimcard,'TOLF',tolf,1.0D-2)
1904 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1905 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1906 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1907 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1908 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1909 & 'Options in energy minimization:'
1910 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1911 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1912 & 'MinMin:',MinMin,' MinFun:',MinFun,
1913 & ' TolF:',TolF,' RTolF:',RTolF
1916 c----------------------------------------------------------------------------
1917 subroutine read_angles(kanal,*)
1919 include 'DIMENSIONS'
1920 include 'COMMON.GEO'
1921 include 'COMMON.VAR'
1922 include 'COMMON.CHAIN'
1923 include 'COMMON.IOUNITS'
1924 include 'COMMON.CONTROL'
1926 c Read angles from input
1928 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1929 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1930 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1931 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1934 c 9/7/01 avoid 180 deg valence angle
1935 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1937 theta(i)=deg2rad*theta(i)
1938 phi(i)=deg2rad*phi(i)
1939 alph(i)=deg2rad*alph(i)
1940 omeg(i)=deg2rad*omeg(i)
1945 c----------------------------------------------------------------------------
1946 subroutine reada(rekord,lancuch,wartosc,default)
1948 character*(*) rekord,lancuch
1949 double precision wartosc,default
1952 iread=index(rekord,lancuch)
1953 if (iread.eq.0) then
1957 iread=iread+ilen(lancuch)+1
1958 read (rekord(iread:),*,err=10,end=10) wartosc
1963 c----------------------------------------------------------------------------
1964 subroutine readi(rekord,lancuch,wartosc,default)
1966 character*(*) rekord,lancuch
1967 integer wartosc,default
1970 iread=index(rekord,lancuch)
1971 if (iread.eq.0) then
1975 iread=iread+ilen(lancuch)+1
1976 read (rekord(iread:),*,err=10,end=10) wartosc
1981 c----------------------------------------------------------------------------
1982 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1985 integer tablica(dim),default
1986 character*(*) rekord,lancuch
1993 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1994 if (iread.eq.0) return
1995 iread=iread+ilen(lancuch)+1
1996 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1999 c----------------------------------------------------------------------------
2000 subroutine multreada(rekord,lancuch,tablica,dim,default)
2003 double precision tablica(dim),default
2004 character*(*) rekord,lancuch
2011 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2012 if (iread.eq.0) return
2013 iread=iread+ilen(lancuch)+1
2014 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2017 c----------------------------------------------------------------------------
2018 subroutine openunits
2020 include 'DIMENSIONS'
2024 character*16 form,nodename
2027 include 'COMMON.SETUP'
2028 include 'COMMON.IOUNITS'
2030 include 'COMMON.CONTROL'
2031 integer lenpre,lenpot,ilen,lentmp,npos
2033 character*3 out1file_text,ucase
2038 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2039 call getenv_loc("PREFIX",prefix)
2041 call getenv_loc("POT",pot)
2042 call getenv_loc("DIRTMP",tmpdir)
2043 call getenv_loc("CURDIR",curdir)
2044 call getenv_loc("OUT1FILE",out1file_text)
2045 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2046 out1file_text=ucase(out1file_text)
2047 if (out1file_text(1:1).eq."Y") then
2050 out1file=fg_rank.gt.0
2055 if (lentmp.gt.0) then
2056 write (*,'(80(1h!))')
2057 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2058 write (*,'(80(1h!))')
2059 write (*,*)"All output files will be on node /tmp directory."
2061 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2062 if (me.eq.king) then
2063 write (*,*) "The master node is ",nodename
2064 else if (fg_rank.eq.0) then
2065 write (*,*) "I am the CG slave node ",nodename
2067 write (*,*) "I am the FG slave node ",nodename
2070 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2071 lenpre = lentmp+lenpre+1
2073 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2074 C Get the names and open the input files
2075 #if defined(WINIFL) || defined(WINPGI)
2076 open(1,file=pref_orig(:ilen(pref_orig))//
2077 & '.inp',status='old',readonly,shared)
2078 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2079 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2080 C Get parameter filenames and open the parameter files.
2081 call getenv_loc('BONDPAR',bondname)
2082 open (ibond,file=bondname,status='old',readonly,shared)
2083 call getenv_loc('THETPAR',thetname)
2084 open (ithep,file=thetname,status='old',readonly,shared)
2085 call getenv_loc('ROTPAR',rotname)
2086 open (irotam,file=rotname,status='old',readonly,shared)
2087 call getenv_loc('TORPAR',torname)
2088 open (itorp,file=torname,status='old',readonly,shared)
2089 call getenv_loc('TORDPAR',tordname)
2090 open (itordp,file=tordname,status='old',readonly,shared)
2091 call getenv_loc('FOURIER',fouriername)
2092 open (ifourier,file=fouriername,status='old',readonly,shared)
2093 call getenv_loc('ELEPAR',elename)
2094 open (ielep,file=elename,status='old',readonly,shared)
2095 call getenv_loc('SIDEPAR',sidename)
2096 open (isidep,file=sidename,status='old',readonly,shared)
2097 call getenv_loc('LIPTRANPAR',liptranname)
2098 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2099 #elif (defined CRAY) || (defined AIX)
2100 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2102 c print *,"Processor",myrank," opened file 1"
2103 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2104 c print *,"Processor",myrank," opened file 9"
2105 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2106 C Get parameter filenames and open the parameter files.
2107 call getenv_loc('BONDPAR',bondname)
2108 open (ibond,file=bondname,status='old',action='read')
2109 c print *,"Processor",myrank," opened file IBOND"
2110 call getenv_loc('THETPAR',thetname)
2111 open (ithep,file=thetname,status='old',action='read')
2113 call getenv_loc('THETPARPDB',thetname_pdb)
2114 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2116 c print *,"Processor",myrank," opened file ITHEP"
2117 call getenv_loc('ROTPAR',rotname)
2118 open (irotam,file=rotname,status='old',action='read')
2120 call getenv_loc('ROTPARPDB',rotname_pdb)
2121 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2123 c print *,"Processor",myrank," opened file IROTAM"
2124 call getenv_loc('TORPAR',torname)
2125 open (itorp,file=torname,status='old',action='read')
2126 c print *,"Processor",myrank," opened file ITORP"
2127 call getenv_loc('TORDPAR',tordname)
2128 open (itordp,file=tordname,status='old',action='read')
2129 c print *,"Processor",myrank," opened file ITORDP"
2130 call getenv_loc('SCCORPAR',sccorname)
2131 open (isccor,file=sccorname,status='old',action='read')
2132 c print *,"Processor",myrank," opened file ISCCOR"
2133 call getenv_loc('FOURIER',fouriername)
2134 open (ifourier,file=fouriername,status='old',action='read')
2135 c print *,"Processor",myrank," opened file IFOURIER"
2136 call getenv_loc('ELEPAR',elename)
2137 open (ielep,file=elename,status='old',action='read')
2138 c print *,"Processor",myrank," opened file IELEP"
2139 call getenv_loc('SIDEPAR',sidename)
2140 open (isidep,file=sidename,status='old',action='read')
2141 call getenv_loc('LIPTRANPAR',liptranname)
2142 open (iliptranpar,file=liptranname,status='old',action='read')
2143 c print *,"Processor",myrank," opened file ISIDEP"
2144 c print *,"Processor",myrank," opened parameter files"
2146 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2147 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2148 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2149 C Get parameter filenames and open the parameter files.
2150 call getenv_loc('BONDPAR',bondname)
2151 open (ibond,file=bondname,status='old')
2152 call getenv_loc('THETPAR',thetname)
2153 open (ithep,file=thetname,status='old')
2155 call getenv_loc('THETPARPDB',thetname_pdb)
2156 open (ithep_pdb,file=thetname_pdb,status='old')
2158 call getenv_loc('ROTPAR',rotname)
2159 open (irotam,file=rotname,status='old')
2161 call getenv_loc('ROTPARPDB',rotname_pdb)
2162 open (irotam_pdb,file=rotname_pdb,status='old')
2164 call getenv_loc('TORPAR',torname)
2165 open (itorp,file=torname,status='old')
2166 call getenv_loc('TORDPAR',tordname)
2167 open (itordp,file=tordname,status='old')
2168 call getenv_loc('SCCORPAR',sccorname)
2169 open (isccor,file=sccorname,status='old')
2170 call getenv_loc('FOURIER',fouriername)
2171 open (ifourier,file=fouriername,status='old')
2172 call getenv_loc('ELEPAR',elename)
2173 open (ielep,file=elename,status='old')
2174 call getenv_loc('SIDEPAR',sidename)
2175 open (isidep,file=sidename,status='old')
2176 call getenv_loc('LIPTRANPAR',liptranname)
2177 open (iliptranpar,file=liptranname,status='old')
2179 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2181 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2182 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2183 C Get parameter filenames and open the parameter files.
2184 call getenv_loc('BONDPAR',bondname)
2185 open (ibond,file=bondname,status='old',readonly)
2186 call getenv_loc('THETPAR',thetname)
2187 open (ithep,file=thetname,status='old',readonly)
2188 call getenv_loc('ROTPAR',rotname)
2189 open (irotam,file=rotname,status='old',readonly)
2190 call getenv_loc('TORPAR',torname)
2191 open (itorp,file=torname,status='old',readonly)
2192 call getenv_loc('TORDPAR',tordname)
2193 open (itordp,file=tordname,status='old',readonly)
2194 call getenv_loc('SCCORPAR',sccorname)
2195 open (isccor,file=sccorname,status='old',readonly)
2197 call getenv_loc('THETPARPDB',thetname_pdb)
2198 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2200 call getenv_loc('FOURIER',fouriername)
2201 open (ifourier,file=fouriername,status='old',readonly)
2202 call getenv_loc('ELEPAR',elename)
2203 open (ielep,file=elename,status='old',readonly)
2204 call getenv_loc('SIDEPAR',sidename)
2205 open (isidep,file=sidename,status='old',readonly)
2206 call getenv_loc('LIPTRANPAR',liptranname)
2207 open (iliptranpar,file=liptranname,status='old',action='read')
2209 call getenv_loc('ROTPARPDB',rotname_pdb)
2210 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2214 call getenv_loc('TUBEPAR',tubename)
2215 #if defined(WINIFL) || defined(WINPGI)
2216 open (itube,file=tubename,status='old',readonly,shared)
2217 #elif (defined CRAY) || (defined AIX)
2218 open (itube,file=tubename,status='old',action='read')
2220 open (itube,file=tubename,status='old')
2222 open (itube,file=tubename,status='old',readonly)
2227 C 8/9/01 In the newest version SCp interaction constants are read from a file
2228 C Use -DOLDSCP to use hard-coded constants instead.
2230 call getenv_loc('SCPPAR',scpname)
2231 #if defined(WINIFL) || defined(WINPGI)
2232 open (iscpp,file=scpname,status='old',readonly,shared)
2233 #elif (defined CRAY) || (defined AIX)
2234 open (iscpp,file=scpname,status='old',action='read')
2236 open (iscpp,file=scpname,status='old')
2238 open (iscpp,file=scpname,status='old',readonly)
2241 call getenv_loc('PATTERN',patname)
2242 #if defined(WINIFL) || defined(WINPGI)
2243 open (icbase,file=patname,status='old',readonly,shared)
2244 #elif (defined CRAY) || (defined AIX)
2245 open (icbase,file=patname,status='old',action='read')
2247 open (icbase,file=patname,status='old')
2249 open (icbase,file=patname,status='old',readonly)
2252 C Open output file only for CG processes
2253 c print *,"Processor",myrank," fg_rank",fg_rank
2254 if (fg_rank.eq.0) then
2256 if (nodes.eq.1) then
2259 npos = dlog10(dfloat(nodes-1))+1
2261 if (npos.lt.3) npos=3
2262 write (liczba,'(i1)') npos
2263 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2265 write (liczba,form) me
2266 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2267 & liczba(:ilen(liczba))
2268 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2270 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2272 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2273 & liczba(:ilen(liczba))//'.mol2'
2274 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2275 & liczba(:ilen(liczba))//'.stat'
2277 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2278 & //liczba(:ilen(liczba))//'.stat')
2279 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2282 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2283 & liczba(:ilen(liczba))//'.const'
2288 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2289 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2290 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2291 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2292 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2294 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2296 rest2name=prefix(:ilen(prefix))//'.rst'
2298 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2301 #if defined(AIX) || defined(PGI) || defined(CRAY)
2302 if (me.eq.king .or. .not. out1file) then
2303 open(iout,file=outname,status='unknown')
2305 open(iout,file="/dev/null",status="unknown")
2309 if (fg_rank.gt.0) then
2310 write (liczba,'(i3.3)') myrank/nfgtasks
2311 write (ll,'(bz,i3.3)') fg_rank
2312 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2318 open(igeom,file=intname,status='unknown',position='append')
2319 open(ipdb,file=pdbname,status='unknown')
2320 open(imol2,file=mol2name,status='unknown')
2321 open(istat,file=statname,status='unknown',position='append')
2323 c1out open(iout,file=outname,status='unknown')
2326 if (me.eq.king .or. .not.out1file)
2327 & open(iout,file=outname,status='unknown')
2329 if (fg_rank.gt.0) then
2330 write (liczba,'(i3.3)') myrank/nfgtasks
2331 write (ll,'(bz,i3.3)') fg_rank
2332 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2337 open(igeom,file=intname,status='unknown',access='append')
2338 open(ipdb,file=pdbname,status='unknown')
2339 open(imol2,file=mol2name,status='unknown')
2340 open(istat,file=statname,status='unknown',access='append')
2342 c1out open(iout,file=outname,status='unknown')
2345 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2346 csa_seed=prefix(:lenpre)//'.CSA.seed'
2347 csa_history=prefix(:lenpre)//'.CSA.history'
2348 csa_bank=prefix(:lenpre)//'.CSA.bank'
2349 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2350 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2351 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2352 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2353 csa_int=prefix(:lenpre)//'.int'
2354 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2355 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2356 csa_in=prefix(:lenpre)//'.CSA.in'
2357 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2360 write (iout,'(80(1h-))')
2361 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2362 write (iout,'(80(1h-))')
2363 write (iout,*) "Input file : ",
2364 & pref_orig(:ilen(pref_orig))//'.inp'
2365 write (iout,*) "Output file : ",
2366 & outname(:ilen(outname))
2368 write (iout,*) "Sidechain potential file : ",
2369 & sidename(:ilen(sidename))
2371 write (iout,*) "SCp potential file : ",
2372 & scpname(:ilen(scpname))
2374 write (iout,*) "Electrostatic potential file : ",
2375 & elename(:ilen(elename))
2376 write (iout,*) "Cumulant coefficient file : ",
2377 & fouriername(:ilen(fouriername))
2378 write (iout,*) "Torsional parameter file : ",
2379 & torname(:ilen(torname))
2380 write (iout,*) "Double torsional parameter file : ",
2381 & tordname(:ilen(tordname))
2382 write (iout,*) "SCCOR parameter file : ",
2383 & sccorname(:ilen(sccorname))
2384 write (iout,*) "Bond & inertia constant file : ",
2385 & bondname(:ilen(bondname))
2386 write (iout,*) "Bending-torsion parameter file : ",
2387 & thetname(:ilen(thetname))
2388 write (iout,*) "Rotamer parameter file : ",
2389 & rotname(:ilen(rotname))
2390 write (iout,*) "Thetpdb parameter file : ",
2391 & thetname_pdb(:ilen(thetname_pdb))
2392 write (iout,*) "Threading database : ",
2393 & patname(:ilen(patname))
2395 &write (iout,*)" DIRTMP : ",
2397 write (iout,'(80(1h-))')
2401 c----------------------------------------------------------------------------
2402 subroutine card_concat(card)
2403 implicit real*8 (a-h,o-z)
2404 include 'DIMENSIONS'
2405 include 'COMMON.IOUNITS'
2407 character*80 karta,ucase
2409 read (inp,'(a)') karta
2412 do while (karta(80:80).eq.'&')
2413 card=card(:ilen(card)+1)//karta(:79)
2414 read (inp,'(a)') karta
2417 card=card(:ilen(card)+1)//karta
2420 c------------------------------------------------------------------------------
2423 include 'DIMENSIONS'
2424 include 'COMMON.CHAIN'
2425 include 'COMMON.IOUNITS'
2426 include 'COMMON.CONTROL'
2428 include 'COMMON.QRESTR'
2430 include 'COMMON.LAGRANGE.5diag'
2432 include 'COMMON.LAGRANGE'
2435 open(irest2,file=rest2name,status='unknown')
2436 read(irest2,*) totT,EK,potE,totE,t_bath
2439 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2442 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2445 read (irest2,*) iset
2450 c------------------------------------------------------------------------------
2451 subroutine read_fragments
2453 include 'DIMENSIONS'
2457 include 'COMMON.SETUP'
2458 include 'COMMON.CHAIN'
2459 include 'COMMON.IOUNITS'
2461 include 'COMMON.QRESTR'
2462 include 'COMMON.CONTROL'
2464 read(inp,*) nset,nfrag,npair,nfrag_back
2465 loc_qlike=(nfrag_back.lt.0)
2466 nfrag_back=iabs(nfrag_back)
2467 if(me.eq.king.or..not.out1file)
2468 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2469 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2471 read(inp,*) mset(iset)
2473 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2475 if(me.eq.king.or..not.out1file)
2476 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2477 & ifrag(2,i,iset), qinfrag(i,iset)
2480 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2482 if(me.eq.king.or..not.out1file)
2483 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2484 & ipair(2,i,iset), qinpair(i,iset)
2488 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2489 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2490 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2491 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2492 if(me.eq.king.or..not.out1file)
2493 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2494 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2495 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2496 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2500 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2501 & wfrag_back(3,i,iset),
2502 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2503 if(me.eq.king.or..not.out1file)
2504 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2505 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2511 C---------------------------------------------------------------------------
2512 subroutine read_afminp
2514 include 'DIMENSIONS'
2518 include 'COMMON.SETUP'
2519 include 'COMMON.CONTROL'
2520 include 'COMMON.CHAIN'
2521 include 'COMMON.IOUNITS'
2522 include 'COMMON.SBRIDGE'
2523 character*320 afmcard
2525 c print *, "wchodze"
2526 call card_concat(afmcard)
2527 call readi(afmcard,"BEG",afmbeg,0)
2528 call readi(afmcard,"END",afmend,0)
2529 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2530 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2531 c print *,'FORCE=' ,forceAFMconst
2532 CCCC NOW PROPERTIES FOR AFM
2535 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2537 distafminit=dsqrt(distafminit)
2538 c print *,'initdist',distafminit
2541 c-------------------------------------------------------------------------------
2542 subroutine read_saxs_constr
2544 include 'DIMENSIONS'
2548 include 'COMMON.SETUP'
2549 include 'COMMON.CONTROL'
2550 include 'COMMON.SAXS'
2551 include 'COMMON.CHAIN'
2552 include 'COMMON.IOUNITS'
2553 include 'COMMON.SBRIDGE'
2554 double precision cm(3),cnorm
2557 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2559 if (saxs_mode.eq.0) then
2560 c SAXS distance distribution
2562 read(inp,*) distsaxs(i),Psaxs(i)
2566 Cnorm = Cnorm + Psaxs(i)
2568 write (iout,*) "Cnorm",Cnorm
2570 Psaxs(i)=Psaxs(i)/Cnorm
2572 write (iout,*) "Normalized distance distribution from SAXS"
2574 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2578 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2580 write (iout,*) "Wsaxs0",Wsaxs0
2584 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2591 cm(j)=cm(j)+Csaxs(j,i)
2599 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2602 write (iout,*) "SAXS sphere coordinates"
2604 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2610 c-------------------------------------------------------------------------------
2611 subroutine read_dist_constr
2613 include 'DIMENSIONS'
2617 include 'COMMON.SETUP'
2618 include 'COMMON.CONTROL'
2619 include 'COMMON.CHAIN'
2620 include 'COMMON.IOUNITS'
2621 include 'COMMON.SBRIDGE'
2622 include 'COMMON.INTERACT'
2623 integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2624 integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2625 double precision wfrag_(100),wpair_(1000)
2626 double precision ddjk,dist,dist_cut,fordepthmax
2627 character*5000 controlcard
2628 logical normalize,next
2630 double precision scal_bfac
2631 double precision xlink(4,0:4) /
2633 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2634 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2635 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2636 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2637 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2638 c print *, "WCHODZE"
2639 write (iout,*) "Calling read_dist_constr"
2640 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2642 restr_on_coord=.false.
2651 call card_concat(controlcard)
2652 next = index(controlcard,"NEXT").gt.0
2653 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2654 write (iout,*) "restr_type",restr_type
2655 call readi(controlcard,"NFRAG",nfrag_,0)
2656 call readi(controlcard,"NFRAG",nfrag_,0)
2657 call readi(controlcard,"NPAIR",npair_,0)
2658 call readi(controlcard,"NDIST",ndist_,0)
2659 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2660 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2661 if (restr_type.eq.10)
2662 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2663 if (restr_type.eq.12)
2664 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2665 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2666 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2667 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2668 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2669 normalize = index(controlcard,"NORMALIZE").gt.0
2670 write (iout,*) "WBOLTZD",wboltzd
2671 write (iout,*) "SCAL_PEAK",scal_peak
2672 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2673 write (iout,*) "IFRAG"
2675 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2677 write (iout,*) "IPAIR"
2679 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2681 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2683 & "Distance restraints as generated from reference structure"
2685 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2686 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2687 & ifrag_(2,i)=nstart_sup+nsup-1
2688 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2690 if (wfrag_(i).eq.0.0d0) cycle
2691 do j=ifrag_(1,i),ifrag_(2,i)-1
2692 do k=j+1,ifrag_(2,i)
2693 c write (iout,*) "j",j," k",k
2695 if (restr_type.eq.1) then
2701 forcon(nhpb)=wfrag_(i)
2702 else if (constr_dist.eq.2) then
2703 if (ddjk.le.dist_cut) then
2709 forcon(nhpb)=wfrag_(i)
2711 else if (restr_type.eq.3) then
2717 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2720 if (.not.out1file .or. me.eq.king)
2721 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2722 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2724 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2725 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2731 if (wpair_(i).eq.0.0d0) cycle
2739 do j=ifrag_(1,ii),ifrag_(2,ii)
2740 do k=ifrag_(1,jj),ifrag_(2,jj)
2742 if (restr_type.eq.1) then
2748 forcon(nhpb)=wpair_(i)
2749 else if (constr_dist.eq.2) then
2750 if (ddjk.le.dist_cut) then
2756 forcon(nhpb)=wpair_(i)
2758 else if (restr_type.eq.3) then
2764 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2767 if (.not.out1file .or. me.eq.king)
2768 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2769 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2771 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2772 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2779 write (iout,*) "Distance restraints as read from input"
2781 if (restr_type.eq.12) then
2782 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2783 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2784 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2785 & fordepth_peak(nhpb_peak+1),npeak
2786 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2787 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2788 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2789 c & fordepth_peak(nhpb_peak+1),npeak
2790 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2791 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2792 nhpb_peak=nhpb_peak+1
2793 irestr_type_peak(nhpb_peak)=12
2794 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2797 if (.not.out1file .or. me.eq.king)
2798 & write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2799 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2800 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2801 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2802 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2804 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2805 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2806 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2807 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2808 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2810 if (ibecarb_peak(nhpb_peak).eq.3) then
2811 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2812 else if (ibecarb_peak(nhpb_peak).eq.2) then
2813 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2814 else if (ibecarb_peak(nhpb_peak).eq.1) then
2815 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2816 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2818 else if (restr_type.eq.11) then
2819 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2820 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2821 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2822 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2824 irestr_type(nhpb)=11
2826 if (.not.out1file .or. me.eq.king)
2827 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2828 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2829 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2831 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2832 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2833 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2835 c if (ibecarb(nhpb).gt.0) then
2836 c ihpb(nhpb)=ihpb(nhpb)+nres
2837 c jhpb(nhpb)=jhpb(nhpb)+nres
2839 if (ibecarb(nhpb).eq.3) then
2840 ihpb(nhpb)=ihpb(nhpb)+nres
2841 else if (ibecarb(nhpb).eq.2) then
2842 ihpb(nhpb)=ihpb(nhpb)+nres
2843 else if (ibecarb(nhpb).eq.1) then
2844 ihpb(nhpb)=ihpb(nhpb)+nres
2845 jhpb(nhpb)=jhpb(nhpb)+nres
2847 else if (restr_type.eq.10) then
2848 c Cross-lonk Markov-like potential
2849 call card_concat(controlcard)
2850 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2851 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2853 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2854 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2855 if (index(controlcard,"ZL").gt.0) then
2857 else if (index(controlcard,"ADH").gt.0) then
2859 else if (index(controlcard,"PDH").gt.0) then
2861 else if (index(controlcard,"DSS").gt.0) then
2866 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2867 & xlink(1,link_type))
2868 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2869 & xlink(2,link_type))
2870 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2871 & xlink(3,link_type))
2872 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2873 & xlink(4,link_type))
2874 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2875 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2876 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2877 if (forcon(nhpb+1).le.0.0d0 .or.
2878 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2880 irestr_type(nhpb)=10
2881 if (ibecarb(nhpb).eq.3) then
2882 jhpb(nhpb)=jhpb(nhpb)+nres
2883 else if (ibecarb(nhpb).eq.2) then
2884 ihpb(nhpb)=ihpb(nhpb)+nres
2885 else if (ibecarb(nhpb).eq.1) then
2886 ihpb(nhpb)=ihpb(nhpb)+nres
2887 jhpb(nhpb)=jhpb(nhpb)+nres
2890 if (.not.out1file .or. me.eq.king)
2891 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2892 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2893 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2896 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2897 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2898 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2903 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2904 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2905 if (forcon(nhpb+1).gt.0.0d0) then
2907 if (dhpb1(nhpb).eq.0.0d0) then
2912 if (ibecarb(nhpb).eq.3) then
2913 jhpb(nhpb)=jhpb(nhpb)+nres
2914 else if (ibecarb(nhpb).eq.2) then
2915 ihpb(nhpb)=ihpb(nhpb)+nres
2916 else if (ibecarb(nhpb).eq.1) then
2917 ihpb(nhpb)=ihpb(nhpb)+nres
2918 jhpb(nhpb)=jhpb(nhpb)+nres
2920 if (dhpb(nhpb).eq.0.0d0)
2921 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2924 if (.not.out1file .or. me.eq.king)
2925 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2926 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2928 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2929 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2932 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2933 C if (forcon(nhpb+1).gt.0.0d0) then
2935 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2938 if (restr_type.eq.4) then
2939 write (iout,*) "The BFAC array"
2941 write (iout,'(i5,f10.5)') i,bfac(i)
2944 if (itype(i).eq.ntyp1) cycle
2946 if (itype(j).eq.ntyp1) cycle
2947 if (itype(i).eq.10) then
2952 if (itype(j).eq.10) then
2962 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2965 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2966 ihpb(nhpb)=i+nres*ii
2967 jhpb(nhpb)=j+nres*jj
2968 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2970 if (.not.out1file .or. me.eq.king) then
2971 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2972 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2973 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2977 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2978 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2979 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2989 if (restr_type.eq.5) then
2990 restr_on_coord=.true.
2992 if (itype(i).eq.ntyp1) cycle
2993 bfac(i)=(scal_bfac/bfac(i))**2
3002 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3003 & fordepthmax=fordepth(i)
3006 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3011 c-------------------------------------------------------------------------------
3012 subroutine read_constr_homology
3014 include 'DIMENSIONS'
3018 include 'COMMON.SETUP'
3019 include 'COMMON.CONTROL'
3020 include 'COMMON.HOMOLOGY'
3021 include 'COMMON.CHAIN'
3022 include 'COMMON.IOUNITS'
3024 include 'COMMON.QRESTR'
3025 include 'COMMON.GEO'
3026 include 'COMMON.INTERACT'
3027 include 'COMMON.NAMES'
3029 c For new homol impl
3031 include 'COMMON.VAR'
3034 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3036 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3037 c & sigma_odl_temp(maxres,maxres,max_template)
3039 character*24 model_ki_dist, model_ki_angle
3040 character*500 controlcard
3041 integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
3042 & ik,iistart,nres_temp
3045 logical liiflag,lfirst
3048 c FP - Nov. 2014 Temporary specifications for new vars
3050 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3052 double precision, dimension (max_template,maxres) :: rescore
3053 double precision, dimension (max_template,maxres) :: rescore2
3054 double precision, dimension (max_template,maxres) :: rescore3
3055 double precision distal
3056 character*24 pdbfile,tpl_k_rescore
3057 c -----------------------------------------------------------------
3058 c Reading multiple PDB ref structures and calculation of retraints
3059 c not using pre-computed ones stored in files model_ki_{dist,angle}
3061 c -----------------------------------------------------------------
3064 c Alternative: reading from input
3065 call card_concat(controlcard)
3066 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3067 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3068 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3069 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3070 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3071 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3072 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3073 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3074 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3075 if(.not.read2sigma.and.start_from_model) then
3076 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3077 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3078 start_from_model=.false.
3079 iranconf=(indpdb.le.0)
3081 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3082 & write(iout,*) 'START_FROM_MODELS is ON'
3083 c if(start_from_model .and. rest) then
3084 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3085 c write(iout,*) 'START_FROM_MODELS is OFF'
3086 c write(iout,*) 'remove restart keyword from input'
3089 if (start_from_model) nmodel_start=constr_homology
3090 if (homol_nset.gt.1)then
3091 call card_concat(controlcard)
3092 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3093 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3094 write(iout,*) "iset homology_weight "
3096 write(iout,*) i,waga_homology(i)
3099 iset=mod(kolor,homol_nset)+1
3102 waga_homology(1)=1.0
3105 cd write (iout,*) "nnt",nnt," nct",nct
3112 c write(iout,*) 'nnt=',nnt,'nct=',nct
3115 do k=1,constr_homology
3128 if (read_homol_frag) then
3129 call read_klapaucjusz
3132 do k=1,constr_homology
3134 read(inp,'(a)') pdbfile
3135 if(me.eq.king .or. .not. out1file)
3136 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3137 & pdbfile(:ilen(pdbfile))
3138 open(ipdbin,file=pdbfile,status='old',err=33)
3140 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3141 & pdbfile(:ilen(pdbfile))
3144 c print *,'Begin reading pdb data'
3146 c Files containing res sim or local scores (former containing sigmas)
3149 write(kic2,'(bz,i2.2)') k
3151 tpl_k_rescore="template"//kic2//".sco"
3155 if (read2sigma) then
3156 call readpdb_template(k)
3163 c Distance restraints
3166 C Copy the coordinates from reference coordinates (?)
3167 do i=1,2*nres_chomo(k)
3170 c write (iout,*) "c(",j,i,") =",c(j,i)
3174 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3176 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3177 open (ientin,file=tpl_k_rescore,status='old')
3178 if (nnt.gt.1) rescore(k,1)=0.0d0
3179 do irec=nnt,nct ! loop for reading res sim
3180 if (read2sigma) then
3181 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3182 & rescore3_tmp,idomain_tmp
3184 idomain(k,i_tmp)=idomain_tmp
3185 rescore(k,i_tmp)=rescore_tmp
3186 rescore2(k,i_tmp)=rescore2_tmp
3187 rescore3(k,i_tmp)=rescore3_tmp
3188 if (.not. out1file .or. me.eq.king)
3189 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3190 & i_tmp,rescore2_tmp,rescore_tmp,
3191 & rescore3_tmp,idomain_tmp
3194 read (ientin,*,end=1401) rescore_tmp
3196 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3197 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3198 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3203 if (waga_dist.ne.0.0d0) then
3211 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3212 c write (iout,*) k,i,j,distal,dist2_cut
3214 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3215 & .and. distal.le.dist2_cut ) then
3221 c write (iout,*) "k",k
3222 c write (iout,*) "i",i," j",j," constr_homology",
3227 if (read2sigma) then
3230 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3232 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3233 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3234 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3236 if (odl(k,ii).le.dist_cut) then
3237 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3240 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3241 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3243 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3244 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3248 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3251 l_homo(k,ii)=.false.
3257 c write (iout,*) "Distance restraints set"
3260 c Theta, dihedral and SC retraints
3262 if (waga_angle.gt.0.0d0) then
3263 c open (ientin,file=tpl_k_sigma_dih,status='old')
3264 c do irec=1,maxres-3 ! loop for reading sigma_dih
3265 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3266 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3267 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3268 c & sigma_dih(k,i+nnt-1)
3273 if (idomain(k,i).eq.0) then
3277 dih(k,i)=phiref(i) ! right?
3278 c read (ientin,*) sigma_dih(k,i) ! original variant
3279 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3280 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3281 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3282 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3283 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3285 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3286 & rescore(k,i-2)+rescore(k,i-3))/4.0
3287 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3288 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3289 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3290 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3291 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3292 c Instead of res sim other local measure of b/b str reliability possible
3293 if (sigma_dih(k,i).ne.0)
3294 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3295 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3299 c write (iout,*) "Dihedral angle restraints set"
3302 if (waga_theta.gt.0.0d0) then
3303 c open (ientin,file=tpl_k_sigma_theta,status='old')
3304 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3305 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3306 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3307 c & sigma_theta(k,i+nnt-1)
3312 do i = nnt+2,nct ! right? without parallel.
3313 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3314 c do i=ithet_start,ithet_end ! with FG parallel.
3315 if (idomain(k,i).eq.0) then
3316 sigma_theta(k,i)=0.0
3319 thetatpl(k,i)=thetaref(i)
3320 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3321 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3322 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3323 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3324 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3325 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3326 & rescore(k,i-2))/3.0
3327 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3328 if (sigma_theta(k,i).ne.0)
3329 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3331 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3332 c rescore(k,i-2) ! right expression ?
3333 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3336 c write (iout,*) "Angle restraints set"
3339 if (waga_d.gt.0.0d0) then
3340 c open (ientin,file=tpl_k_sigma_d,status='old')
3341 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3342 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3343 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3344 c & sigma_d(k,i+nnt-1)
3348 do i = nnt,nct ! right? without parallel.
3349 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3350 c do i=loc_start,loc_end ! with FG parallel.
3351 if (itype(i).eq.10) cycle
3352 if (idomain(k,i).eq.0 ) then
3359 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3360 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3361 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3362 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3363 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3364 if (sigma_d(k,i).ne.0)
3365 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3367 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3368 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3369 c read (ientin,*) sigma_d(k,i) ! 1st variant
3373 c write (iout,*) "SC restraints set"
3376 c remove distance restraints not used in any model from the list
3377 c shift data in all arrays
3379 c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
3380 if (waga_dist.ne.0.0d0) then
3387 c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3388 c & .and. distal.le.dist2_cut ) then
3389 c write (iout,*) "i",i," j",j," ii",ii
3391 if (ii_in_use(ii).eq.0.and.liiflag.or.
3392 & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3399 if(i10.eq.lim_odl) i10=i10+1
3401 ires_homo(iistart+ki)=ires_homo(ki+i01)
3402 jres_homo(iistart+ki)=jres_homo(ki+i01)
3403 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3404 do k=1,constr_homology
3405 odl(k,iistart+ki)=odl(k,ki+i01)
3406 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3407 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3410 iistart=iistart+i10-i01
3413 if (ii_in_use(ii).ne.0.and..not.liiflag) then
3421 c write (iout,*) "Removing distances completed"
3423 endif ! .not. klapaucjusz
3425 if (constr_homology.gt.0) call homology_partition
3426 c write (iout,*) "After homology_partition"
3428 if (constr_homology.gt.0) call init_int_table
3429 c write (iout,*) "After init_int_table"
3431 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3432 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3436 if (.not.out_template_restr) return
3437 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3438 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3439 write (iout,*) "Distance restraints from templates"
3441 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3442 & ii,ires_homo(ii),jres_homo(ii),
3443 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3444 & ki=1,constr_homology)
3446 write (iout,*) "Dihedral angle restraints from templates"
3448 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3449 & (rad2deg*dih(ki,i),
3450 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3452 write (iout,*) "Virtual-bond angle restraints from templates"
3454 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3455 & (rad2deg*thetatpl(ki,i),
3456 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3458 write (iout,*) "SC restraints from templates"
3460 write(iout,'(i5,100(4f8.2,4x))') i,
3461 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3462 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3465 c -----------------------------------------------------------------
3468 c----------------------------------------------------------------------
3470 subroutine flush(iu)
3475 subroutine flush(iu)
3480 c------------------------------------------------------------------------------
3481 subroutine copy_to_tmp(source)
3483 include "DIMENSIONS"
3484 include "COMMON.IOUNITS"
3485 character*(*) source
3486 character* 256 tmpfile
3490 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3491 inquire(file=tmpfile,exist=ex)
3493 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3494 & " to temporary directory..."
3495 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3496 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3500 c------------------------------------------------------------------------------
3501 subroutine move_from_tmp(source)
3503 include "DIMENSIONS"
3504 include "COMMON.IOUNITS"
3505 character*(*) source
3508 write (*,*) "Moving ",source(:ilen(source)),
3509 & " from temporary directory to working directory"
3510 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3511 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3514 c------------------------------------------------------------------------------
3515 subroutine random_init(seed)
3517 C Initialize random number generator
3520 include 'DIMENSIONS'
3523 logical OKRandom, prng_restart
3525 integer iseed_array(4)
3526 integer error_msg,ierr
3528 include 'COMMON.IOUNITS'
3529 include 'COMMON.TIME1'
3530 include 'COMMON.THREAD'
3531 include 'COMMON.SBRIDGE'
3532 include 'COMMON.CONTROL'
3533 include 'COMMON.MCM'
3534 include 'COMMON.MAP'
3535 include 'COMMON.HEADER'
3536 include 'COMMON.CSA'
3537 include 'COMMON.CHAIN'
3538 include 'COMMON.MUCA'
3540 include 'COMMON.FFIELD'
3541 include 'COMMON.SETUP'
3543 double precision seed,ran_number
3544 iseed=-dint(dabs(seed))
3545 if (iseed.eq.0) then
3546 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3547 & 'Random seed undefined. The program will stop.'
3548 write (*,'(/80(1h*)/20x,a/80(1h*))')
3549 & 'Random seed undefined. The program will stop.'
3551 call mpi_finalize(mpi_comm_world,ierr)
3553 stop 'Bad random seed.'
3556 if (fg_rank.eq.0) then
3560 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3561 OKRandom = prng_restart(me,iseed)
3564 tmp=65536.0d0**(4-i)
3565 iseed_array(i) = dint(seed/tmp)
3566 seed=seed-iseed_array(i)*tmp
3569 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3570 & (iseed_array(i),i=1,4)
3571 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3572 & (iseed_array(i),i=1,4)
3573 OKRandom = prng_restart(me,iseed_array)
3576 c r1 = prng_next(me)
3577 r1=ran_number(0.0D0,1.0D0)
3579 & write (iout,*) 'ran_num',r1
3580 if (r1.lt.0.0d0) OKRandom=.false.
3582 if (.not.OKRandom) then
3583 write (iout,*) 'PRNG IS NOT WORKING!!!'
3584 print *,'PRNG IS NOT WORKING!!!'
3587 call mpi_abort(mpi_comm_world,error_msg,ierr)
3590 write (iout,*) 'too many processors for parallel prng'
3591 write (*,*) 'too many processors for parallel prng'
3599 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3603 c----------------------------------------------------------------------
3604 subroutine read_klapaucjusz
3606 include 'DIMENSIONS'
3610 include 'COMMON.SETUP'
3611 include 'COMMON.CONTROL'
3612 include 'COMMON.HOMOLOGY'
3613 include 'COMMON.CHAIN'
3614 include 'COMMON.IOUNITS'
3616 include 'COMMON.GEO'
3617 include 'COMMON.INTERACT'
3618 include 'COMMON.NAMES'
3619 character*256 fragfile
3620 integer ninclust(maxclust),inclust(max_template,maxclust),
3621 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
3624 character*24 model_ki_dist, model_ki_angle
3625 character*500 controlcard
3626 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3627 & ik,ll,ii,kk,iistart,iishift,lim_xx
3628 double precision distal
3629 logical lprn /.true./
3636 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3637 double precision, dimension (max_template,maxres) :: rescore
3638 double precision, dimension (max_template,maxres) :: rescore2
3639 character*24 pdbfile,tpl_k_rescore
3642 c For new homol impl
3644 include 'COMMON.VAR'
3646 call getenv("FRAGFILE",fragfile)
3647 open(ientin,file=fragfile,status="old",err=10)
3648 read(ientin,*) constr_homology,nclust
3654 do k=1,constr_homology
3655 read(ientin,'(a)') pdbfile
3656 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3657 & pdbfile(:ilen(pdbfile))
3658 open(ipdbin,file=pdbfile,status='old',err=33)
3660 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3661 & pdbfile(:ilen(pdbfile))
3666 call readpdb_template(k)
3676 read(ientin,*) ninclust(i),nresclust(i)
3677 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3678 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3681 c Loop over clusters
3684 do ll = 1,ninclust(l)
3692 idomain(k,iresclust(i,l)+1) = 1
3694 idomain(k,iresclust(i,l)) = 1
3698 c Distance restraints
3701 C Copy the coordinates from reference coordinates (?)
3707 c write (iout,*) "c(",j,i,") =",c(j,i)
3710 call int_from_cart(.true.,.false.)
3711 call sc_loc_geom(.false.)
3713 thetaref(i)=theta(i)
3717 if (waga_dist.ne.0.0d0) then
3725 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3726 c write (iout,*) k,i,j,distal,dist2_cut
3728 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3729 & .and. distal.le.dist2_cut ) then
3735 c write (iout,*) "k",k
3736 c write (iout,*) "i",i," j",j," constr_homology",
3741 if (read2sigma) then
3744 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3746 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3747 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3748 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3750 if (odl(k,ii).le.dist_cut) then
3751 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3754 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3755 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3757 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3758 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3762 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3765 c l_homo(k,ii)=.false.
3772 c Theta, dihedral and SC retraints
3774 if (waga_angle.gt.0.0d0) then
3776 if (idomain(k,i).eq.0) then
3777 c sigma_dih(k,i)=0.0
3781 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3782 & rescore(k,i-2)+rescore(k,i-3))/4.0
3783 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3784 c & " sigma_dihed",sigma_dih(k,i)
3785 if (sigma_dih(k,i).ne.0)
3786 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3791 if (waga_theta.gt.0.0d0) then
3793 if (idomain(k,i).eq.0) then
3794 c sigma_theta(k,i)=0.0
3797 thetatpl(k,i)=thetaref(i)
3798 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3799 & rescore(k,i-2))/3.0
3800 if (sigma_theta(k,i).ne.0)
3801 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3805 if (waga_d.gt.0.0d0) then
3807 if (itype(i).eq.10) cycle
3808 if (idomain(k,i).eq.0 ) then
3815 sigma_d(k,i)=rescore(k,i)
3816 if (sigma_d(k,i).ne.0)
3817 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3818 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3824 c remove distance restraints not used in any model from the list
3825 c shift data in all arrays
3827 if (waga_dist.ne.0.0d0) then
3833 if (ii_in_use(ii).eq.0.and.liiflag) then
3837 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3838 & .not.liiflag.and.ii.eq.lim_odl) then
3839 if (ii.eq.lim_odl) then
3840 iishift=ii-iistart+1
3845 do ki=iistart,lim_odl-iishift
3846 ires_homo(ki)=ires_homo(ki+iishift)
3847 jres_homo(ki)=jres_homo(ki+iishift)
3848 ii_in_use(ki)=ii_in_use(ki+iishift)
3849 do k=1,constr_homology
3850 odl(k,ki)=odl(k,ki+iishift)
3851 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3852 l_homo(k,ki)=l_homo(k,ki+iishift)
3856 lim_odl=lim_odl-iishift
3863 10 stop "Error in fragment file"