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))
1188 call readpdb_template(nmodel_start+1)
1190 if (nres.ge.nres_temp) then
1191 nmodel_start=nmodel_start+1
1192 pdbfiles_chomo(nmodel_start)=pdbfile
1195 chomo(j,i,nmodel_start)=c(j,i)
1199 if (me.eq.king .or. .not. out1file)
1200 & write (iout,'(a,2i5,1x,a)')
1201 & "Different number of residues",nres_temp,nres,
1207 if (nmodel_start.eq.0) then
1208 if (me.eq.king .or. .not. out1file)
1209 & write (iout,'(a)')
1210 & "No valid starting model found START_FROM_MODELS is OFF"
1211 start_from_model=.false.
1213 write (iout,*) "nmodel_start",nmodel_start
1219 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1220 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1221 & modecalc.ne.10) then
1222 C If input structure hasn't been supplied from the PDB file read or generate
1224 if (iranconf.eq.0 .and. .not. extconf .and. .not.
1225 & start_from_model) then
1226 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1227 & write (iout,'(a)') 'Initial geometry will be read in.'
1229 read(inp,'(8f10.5)',end=36,err=36)
1230 & ((c(l,k),l=1,3),k=1,nres),
1231 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1232 if (nnt.gt.1) c(:,nres+1)=c(:,1)
1233 if (nct.lt.nres) c(:,2*nres)=c(:,nres)
1234 c write (iout,*) "Exit READ_CART"
1235 c write (iout,'(8f10.5)')
1236 c & ((c(l,k),l=1,3),k=1,nres),
1237 c & ((c(l,k+nres),l=1,3),k=nnt,nct)
1244 dc(j,i)=c(j,i+1)-c(j,i)
1245 c dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1249 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1251 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1252 c dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1257 c dc_norm(j,i+nres)=0.0d0
1261 call int_from_cart1(.true.)
1262 write (iout,*) "Finish INT_TO_CART"
1263 c write (iout,*) "DC"
1265 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1266 c & (dc(j,i+nres),j=1,3)
1272 call read_angles(inp,*36)
1274 call chainbuild_extconf
1277 36 write (iout,'(a)') 'Error reading angle file.'
1279 call mpi_finalize( MPI_COMM_WORLD,IERR )
1281 stop 'Error reading angle file.'
1283 else if (extconf) then
1284 if (me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1285 & write (iout,'(a)') 'Extended chain initial geometry.'
1287 theta(i)=90d0*deg2rad
1290 phi(i)=180d0*deg2rad
1293 alph(i)=110d0*deg2rad
1296 omeg(i)=-120d0*deg2rad
1297 if (itype(i).le.0) omeg(i)=-omeg(i)
1300 call chainbuild_extconf
1301 else if (.not. start_from_model) then
1302 if(me.eq.king.or..not.out1file)
1303 & write (iout,'(a)') 'Random-generated initial geometry.'
1306 if (me.eq.king .or. fg_rank.eq.0 .and. (
1307 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1311 call gen_rand_conf(itmp,*30)
1313 30 write (iout,*) 'Failed to generate random conformation',
1314 & ', itrial=',itrial
1315 write (*,*) 'Processor:',me,
1316 & ' Failed to generate random conformation',
1326 write (iout,'(a,i3,a)') 'Processor:',me,
1327 & ' error in generating random conformation.'
1328 write (*,'(a,i3,a)') 'Processor:',me,
1329 & ' error in generating random conformation.'
1332 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1337 & ' error in generating random conformation.'
1342 elseif (modecalc.eq.4) then
1343 read (inp,'(a)') intinname
1344 open (intin,file=intinname,status='old',err=333)
1345 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1346 & write (iout,'(a)') 'intinname',intinname
1347 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1349 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1351 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1353 stop 'Error opening angle file.'
1357 C Generate distance constraints, if the PDB structure is to be regularized.
1358 if (nthread.gt.0) then
1359 call read_threadbase
1362 if (me.eq.king .or. .not. out1file)
1364 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1365 write (iout,'(/a,i3,a)')
1366 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1367 write (iout,'(20i4)') (iss(i),i=1,ns)
1369 write(iout,*)"Running with dynamic disulfide-bond formation"
1371 write (iout,'(/a/)') 'Pre-formed links are:'
1377 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1378 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1384 if (ns.gt.0.and.dyn_ss) then
1388 forcon(i-nss)=forcon(i)
1395 dyn_ss_mask(iss(i))=.true.
1400 c write (iout,*) "DC"
1402 c write (iout,'(i5,3f10.5,5x,3f10.5)') i,(dc(j,i),j=1,3),
1403 c & (dc(j,i+nres),j=1,3)
1405 if (i2ndstr.gt.0) call secstrp2dihc
1406 c call geom_to_var(nvar,x)
1407 c call etotal(energia(0))
1408 c call enerprint(energia(0))
1409 c call briefout(0,etot)
1411 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1412 cd write (iout,'(a)') 'Variable list:'
1413 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1415 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1416 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1417 & 'Processor',myrank,': end reading molecular data.'
1422 c--------------------------------------------------------------------------
1423 logical function seq_comp(itypea,itypeb,length)
1425 integer length,itypea(length),itypeb(length)
1428 if (itypea(i).ne.itypeb(i)) then
1436 c-----------------------------------------------------------------------------
1437 subroutine read_bridge
1438 C Read information about disulfide bridges.
1440 include 'DIMENSIONS'
1445 include 'COMMON.IOUNITS'
1446 include 'COMMON.GEO'
1447 include 'COMMON.VAR'
1448 include 'COMMON.INTERACT'
1449 include 'COMMON.LOCAL'
1450 include 'COMMON.NAMES'
1451 include 'COMMON.CHAIN'
1452 include 'COMMON.FFIELD'
1453 include 'COMMON.SBRIDGE'
1454 include 'COMMON.HEADER'
1455 include 'COMMON.CONTROL'
1456 include 'COMMON.DBASE'
1457 include 'COMMON.THREAD'
1458 include 'COMMON.TIME1'
1459 include 'COMMON.SETUP'
1461 C Read bridging residues.
1462 read (inp,*) ns,(iss(i),i=1,ns)
1463 c 5/24/2020 Adam: Added a table to translate residue numbers to cysteine
1470 if(me.eq.king.or..not.out1file)
1471 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1472 C Check whether the specified bridging residues are cystines.
1474 if (itype(iss(i)).ne.1) then
1475 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1476 & 'Do you REALLY think that the residue ',
1477 & restyp(itype(iss(i))),i,
1478 & ' can form a disulfide bridge?!!!'
1479 write (*,'(2a,i3,a)')
1480 & 'Do you REALLY think that the residue ',
1481 & restyp(itype(iss(i))),i,
1482 & ' can form a disulfide bridge?!!!'
1484 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1489 C Read preformed bridges.
1491 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1493 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1496 C Check if the residues involved in bridges are in the specified list of
1497 C bridging residues.
1500 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1501 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1502 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1503 & ' contains residues present in other pairs.'
1504 write (*,'(a,i3,a)') 'Disulfide pair',i,
1505 & ' contains residues present in other pairs.'
1507 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1513 if (ihpb(i).eq.iss(j)) goto 10
1515 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1518 if (jhpb(i).eq.iss(j)) goto 20
1520 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1526 ihpb(i)=ihpb(i)+nres
1527 jhpb(i)=jhpb(i)+nres
1533 c----------------------------------------------------------------------------
1534 subroutine read_x(kanal,*)
1536 include 'DIMENSIONS'
1537 include 'COMMON.GEO'
1538 include 'COMMON.VAR'
1539 include 'COMMON.CHAIN'
1540 include 'COMMON.IOUNITS'
1541 include 'COMMON.CONTROL'
1542 include 'COMMON.LOCAL'
1543 include 'COMMON.INTERACT'
1544 integer i,j,k,l,kanal
1545 c Read coordinates from input
1547 read(kanal,'(8f10.5)',end=10,err=10)
1548 & ((c(l,k),l=1,3),k=1,nres),
1549 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1552 c(j,2*nres)=c(j,nres)
1554 call int_from_cart1(.false.)
1557 dc(j,i)=c(j,i+1)-c(j,i)
1558 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1562 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1564 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1565 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1573 c----------------------------------------------------------------------------
1574 subroutine read_threadbase
1576 include 'DIMENSIONS'
1577 include 'COMMON.IOUNITS'
1578 include 'COMMON.GEO'
1579 include 'COMMON.VAR'
1580 include 'COMMON.INTERACT'
1581 include 'COMMON.LOCAL'
1582 include 'COMMON.NAMES'
1583 include 'COMMON.CHAIN'
1584 include 'COMMON.FFIELD'
1585 include 'COMMON.SBRIDGE'
1586 include 'COMMON.HEADER'
1587 include 'COMMON.CONTROL'
1588 include 'COMMON.DBASE'
1589 include 'COMMON.THREAD'
1590 include 'COMMON.TIME1'
1592 double precision dist
1593 C Read pattern database for threading.
1594 read (icbase,*) nseq
1596 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1597 & nres_base(2,i),nres_base(3,i)
1598 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1600 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1601 c & nres_base(2,i),nres_base(3,i)
1602 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1606 if (weidis.eq.0.0D0) weidis=0.1D0
1615 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1616 write (iout,'(a,i5)') 'nexcl: ',nexcl
1617 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1620 c------------------------------------------------------------------------------
1621 subroutine setup_var
1623 include 'DIMENSIONS'
1624 include 'COMMON.IOUNITS'
1625 include 'COMMON.GEO'
1626 include 'COMMON.VAR'
1627 include 'COMMON.INTERACT'
1628 include 'COMMON.LOCAL'
1629 include 'COMMON.NAMES'
1630 include 'COMMON.CHAIN'
1631 include 'COMMON.FFIELD'
1632 include 'COMMON.SBRIDGE'
1633 include 'COMMON.HEADER'
1634 include 'COMMON.CONTROL'
1635 include 'COMMON.DBASE'
1636 include 'COMMON.THREAD'
1637 include 'COMMON.TIME1'
1639 C Set up variable list.
1645 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1647 ialph(i,1)=nvar+nside
1651 if (indphi.gt.0) then
1653 else if (indback.gt.0) then
1660 c----------------------------------------------------------------------------
1661 subroutine gen_dist_constr
1662 C Generate CA distance constraints.
1664 include 'DIMENSIONS'
1665 include 'COMMON.IOUNITS'
1666 include 'COMMON.GEO'
1667 include 'COMMON.VAR'
1668 include 'COMMON.INTERACT'
1669 include 'COMMON.LOCAL'
1670 include 'COMMON.NAMES'
1671 include 'COMMON.CHAIN'
1672 include 'COMMON.FFIELD'
1673 include 'COMMON.SBRIDGE'
1674 include 'COMMON.HEADER'
1675 include 'COMMON.CONTROL'
1676 include 'COMMON.DBASE'
1677 include 'COMMON.THREAD'
1678 include 'COMMON.SPLITELE'
1679 include 'COMMON.TIME1'
1680 integer i,j,itype_pdb(maxres)
1681 common /pizda/ itype_pdb
1683 double precision dist
1685 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1686 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1687 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1689 do i=nstart_sup,nstart_sup+nsup-1
1690 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1691 cd & ' seq_pdb', restyp(itype_pdb(i))
1692 do j=i+2,nstart_sup+nsup-1
1693 c 5/24/2020 Adam: Cutoff included to reduce array size
1695 if (dd.gt.r_cut_int) cycle
1697 ihpb(nhpb)=i+nstart_seq-nstart_sup
1698 jhpb(nhpb)=j+nstart_seq-nstart_sup
1703 cd write (iout,'(a)') 'Distance constraints:'
1708 cd if (ii.gt.nres) then
1713 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1714 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1715 cd & dhpb(i),forcon(i)
1719 c----------------------------------------------------------------------------
1722 include 'DIMENSIONS'
1723 include 'COMMON.MAP'
1724 include 'COMMON.IOUNITS'
1726 character*3 angid(4) /'THE','PHI','ALP','OME'/
1727 character*80 mapcard,ucase
1729 read (inp,'(a)') mapcard
1730 mapcard=ucase(mapcard)
1731 if (index(mapcard,'PHI').gt.0) then
1733 else if (index(mapcard,'THE').gt.0) then
1735 else if (index(mapcard,'ALP').gt.0) then
1737 else if (index(mapcard,'OME').gt.0) then
1740 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1741 stop 'Error - illegal variable spec in MAP card.'
1743 call readi (mapcard,'RES1',res1(imap),0)
1744 call readi (mapcard,'RES2',res2(imap),0)
1745 if (res1(imap).eq.0) then
1746 res1(imap)=res2(imap)
1747 else if (res2(imap).eq.0) then
1748 res2(imap)=res1(imap)
1750 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1752 & 'Error - illegal definition of variable group in MAP.'
1753 stop 'Error - illegal definition of variable group in MAP.'
1755 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1756 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1757 call readi(mapcard,'NSTEP',nstep(imap),0)
1758 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1760 & 'Illegal boundary and/or step size specification in MAP.'
1761 stop 'Illegal boundary and/or step size specification in MAP.'
1766 c----------------------------------------------------------------------------
1769 include 'DIMENSIONS'
1770 include 'COMMON.IOUNITS'
1771 include 'COMMON.GEO'
1772 include 'COMMON.CSA'
1773 include 'COMMON.BANK'
1774 include 'COMMON.CONTROL'
1776 character*620 mcmcard
1777 call card_concat(mcmcard)
1779 call readi(mcmcard,'NCONF',nconf,50)
1780 call readi(mcmcard,'NADD',nadd,0)
1781 call readi(mcmcard,'JSTART',jstart,1)
1782 call readi(mcmcard,'JEND',jend,1)
1783 call readi(mcmcard,'NSTMAX',nstmax,500000)
1784 call readi(mcmcard,'N0',n0,1)
1785 call readi(mcmcard,'N1',n1,6)
1786 call readi(mcmcard,'N2',n2,4)
1787 call readi(mcmcard,'N3',n3,0)
1788 call readi(mcmcard,'N4',n4,0)
1789 call readi(mcmcard,'N5',n5,0)
1790 call readi(mcmcard,'N6',n6,10)
1791 call readi(mcmcard,'N7',n7,0)
1792 call readi(mcmcard,'N8',n8,0)
1793 call readi(mcmcard,'N9',n9,0)
1794 call readi(mcmcard,'N14',n14,0)
1795 call readi(mcmcard,'N15',n15,0)
1796 call readi(mcmcard,'N16',n16,0)
1797 call readi(mcmcard,'N17',n17,0)
1798 call readi(mcmcard,'N18',n18,0)
1800 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1802 call readi(mcmcard,'NDIFF',ndiff,2)
1803 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1804 call readi(mcmcard,'IS1',is1,1)
1805 call readi(mcmcard,'IS2',is2,8)
1806 call readi(mcmcard,'NRAN0',nran0,4)
1807 call readi(mcmcard,'NRAN1',nran1,2)
1808 call readi(mcmcard,'IRR',irr,1)
1809 call readi(mcmcard,'NSEED',nseed,20)
1810 call readi(mcmcard,'NTOTAL',ntotal,10000)
1811 call reada(mcmcard,'CUT1',cut1,2.0d0)
1812 call reada(mcmcard,'CUT2',cut2,5.0d0)
1813 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1814 call readi(mcmcard,'ICMAX',icmax,3)
1815 call readi(mcmcard,'IRESTART',irestart,0)
1816 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1819 call reada(mcmcard,'DELE',dele,20.0d0)
1820 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1821 call readi(mcmcard,'IREF',iref,0)
1822 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1823 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1824 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1825 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1826 write (iout,*) "NCONF_IN",nconf_in
1829 c----------------------------------------------------------------------------
1832 include 'DIMENSIONS'
1833 include 'COMMON.MCM'
1834 include 'COMMON.MCE'
1835 include 'COMMON.IOUNITS'
1837 character*320 mcmcard
1839 call card_concat(mcmcard)
1840 call readi(mcmcard,'MAXACC',maxacc,100)
1841 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1842 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1843 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1844 call readi(mcmcard,'MAXREPM',maxrepm,200)
1845 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1846 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1847 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1848 call reada(mcmcard,'E_UP',e_up,5.0D0)
1849 call reada(mcmcard,'DELTE',delte,0.1D0)
1850 call readi(mcmcard,'NSWEEP',nsweep,5)
1851 call readi(mcmcard,'NSTEPH',nsteph,0)
1852 call readi(mcmcard,'NSTEPC',nstepc,0)
1853 call reada(mcmcard,'TMIN',tmin,298.0D0)
1854 call reada(mcmcard,'TMAX',tmax,298.0D0)
1855 call readi(mcmcard,'NWINDOW',nwindow,0)
1856 call readi(mcmcard,'PRINT_MC',print_mc,0)
1857 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1858 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1859 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1860 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1861 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1862 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1863 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1864 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1865 if (nwindow.gt.0) then
1866 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1868 winlen(i)=winend(i)-winstart(i)+1
1871 if (tmax.lt.tmin) tmax=tmin
1872 if (tmax.eq.tmin) then
1876 if (nstepc.gt.0 .and. nsteph.gt.0) then
1877 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1878 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1880 C Probabilities of different move types
1881 sumpro_type(0)=0.0D0
1882 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1883 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1884 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1885 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1886 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1887 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1888 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1890 print *,'i',i,' sumprotype',sumpro_type(i)
1891 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1892 print *,'i',i,' sumprotype',sumpro_type(i)
1896 c----------------------------------------------------------------------------
1897 subroutine read_minim
1899 include 'DIMENSIONS'
1900 include 'COMMON.MINIM'
1901 include 'COMMON.IOUNITS'
1903 character*320 minimcard
1904 call card_concat(minimcard)
1905 call readi(minimcard,'MAXMIN',maxmin,2000)
1906 call readi(minimcard,'MAXFUN',maxfun,5000)
1907 call readi(minimcard,'MINMIN',minmin,maxmin)
1908 call readi(minimcard,'MINFUN',minfun,maxmin)
1909 call reada(minimcard,'TOLF',tolf,1.0D-2)
1910 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1911 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1912 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1913 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1914 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1915 & 'Options in energy minimization:'
1916 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1917 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1918 & 'MinMin:',MinMin,' MinFun:',MinFun,
1919 & ' TolF:',TolF,' RTolF:',RTolF
1922 c----------------------------------------------------------------------------
1923 subroutine read_angles(kanal,*)
1925 include 'DIMENSIONS'
1926 include 'COMMON.GEO'
1927 include 'COMMON.VAR'
1928 include 'COMMON.CHAIN'
1929 include 'COMMON.IOUNITS'
1930 include 'COMMON.CONTROL'
1932 c Read angles from input
1934 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1935 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1936 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1937 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1940 c 9/7/01 avoid 180 deg valence angle
1941 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1943 theta(i)=deg2rad*theta(i)
1944 phi(i)=deg2rad*phi(i)
1945 alph(i)=deg2rad*alph(i)
1946 omeg(i)=deg2rad*omeg(i)
1951 c----------------------------------------------------------------------------
1952 subroutine reada(rekord,lancuch,wartosc,default)
1954 character*(*) rekord,lancuch
1955 double precision wartosc,default
1958 iread=index(rekord,lancuch)
1959 if (iread.eq.0) then
1963 iread=iread+ilen(lancuch)+1
1964 read (rekord(iread:),*,err=10,end=10) wartosc
1969 c----------------------------------------------------------------------------
1970 subroutine readi(rekord,lancuch,wartosc,default)
1972 character*(*) rekord,lancuch
1973 integer wartosc,default
1976 iread=index(rekord,lancuch)
1977 if (iread.eq.0) then
1981 iread=iread+ilen(lancuch)+1
1982 read (rekord(iread:),*,err=10,end=10) wartosc
1987 c----------------------------------------------------------------------------
1988 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1991 integer tablica(dim),default
1992 character*(*) rekord,lancuch
1999 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2000 if (iread.eq.0) return
2001 iread=iread+ilen(lancuch)+1
2002 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2005 c----------------------------------------------------------------------------
2006 subroutine multreada(rekord,lancuch,tablica,dim,default)
2009 double precision tablica(dim),default
2010 character*(*) rekord,lancuch
2017 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2018 if (iread.eq.0) return
2019 iread=iread+ilen(lancuch)+1
2020 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2023 c----------------------------------------------------------------------------
2024 subroutine openunits
2026 include 'DIMENSIONS'
2030 character*16 form,nodename
2033 include 'COMMON.SETUP'
2034 include 'COMMON.IOUNITS'
2036 include 'COMMON.CONTROL'
2037 integer lenpre,lenpot,ilen,lentmp,npos
2039 character*3 out1file_text,ucase
2044 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2045 call getenv_loc("PREFIX",prefix)
2047 call getenv_loc("POT",pot)
2048 call getenv_loc("DIRTMP",tmpdir)
2049 call getenv_loc("CURDIR",curdir)
2050 call getenv_loc("OUT1FILE",out1file_text)
2051 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2052 out1file_text=ucase(out1file_text)
2053 if (out1file_text(1:1).eq."Y") then
2056 out1file=fg_rank.gt.0
2061 if (lentmp.gt.0) then
2062 write (*,'(80(1h!))')
2063 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2064 write (*,'(80(1h!))')
2065 write (*,*)"All output files will be on node /tmp directory."
2067 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2068 if (me.eq.king) then
2069 write (*,*) "The master node is ",nodename
2070 else if (fg_rank.eq.0) then
2071 write (*,*) "I am the CG slave node ",nodename
2073 write (*,*) "I am the FG slave node ",nodename
2076 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2077 lenpre = lentmp+lenpre+1
2079 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2080 C Get the names and open the input files
2081 #if defined(WINIFL) || defined(WINPGI)
2082 open(1,file=pref_orig(:ilen(pref_orig))//
2083 & '.inp',status='old',readonly,shared)
2084 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2085 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2086 C Get parameter filenames and open the parameter files.
2087 call getenv_loc('BONDPAR',bondname)
2088 open (ibond,file=bondname,status='old',readonly,shared)
2089 call getenv_loc('THETPAR',thetname)
2090 open (ithep,file=thetname,status='old',readonly,shared)
2091 call getenv_loc('ROTPAR',rotname)
2092 open (irotam,file=rotname,status='old',readonly,shared)
2093 call getenv_loc('TORPAR',torname)
2094 open (itorp,file=torname,status='old',readonly,shared)
2095 call getenv_loc('TORDPAR',tordname)
2096 open (itordp,file=tordname,status='old',readonly,shared)
2097 call getenv_loc('FOURIER',fouriername)
2098 open (ifourier,file=fouriername,status='old',readonly,shared)
2099 call getenv_loc('ELEPAR',elename)
2100 open (ielep,file=elename,status='old',readonly,shared)
2101 call getenv_loc('SIDEPAR',sidename)
2102 open (isidep,file=sidename,status='old',readonly,shared)
2103 call getenv_loc('LIPTRANPAR',liptranname)
2104 open (iliptranpar,file=liptranname,status='old',readonly,shared)
2105 #elif (defined CRAY) || (defined AIX)
2106 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2108 c print *,"Processor",myrank," opened file 1"
2109 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2110 c print *,"Processor",myrank," opened file 9"
2111 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2112 C Get parameter filenames and open the parameter files.
2113 call getenv_loc('BONDPAR',bondname)
2114 open (ibond,file=bondname,status='old',action='read')
2115 c print *,"Processor",myrank," opened file IBOND"
2116 call getenv_loc('THETPAR',thetname)
2117 open (ithep,file=thetname,status='old',action='read')
2119 call getenv_loc('THETPARPDB',thetname_pdb)
2120 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2122 c print *,"Processor",myrank," opened file ITHEP"
2123 call getenv_loc('ROTPAR',rotname)
2124 open (irotam,file=rotname,status='old',action='read')
2126 call getenv_loc('ROTPARPDB',rotname_pdb)
2127 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2129 c print *,"Processor",myrank," opened file IROTAM"
2130 call getenv_loc('TORPAR',torname)
2131 open (itorp,file=torname,status='old',action='read')
2132 c print *,"Processor",myrank," opened file ITORP"
2133 call getenv_loc('TORDPAR',tordname)
2134 open (itordp,file=tordname,status='old',action='read')
2135 c print *,"Processor",myrank," opened file ITORDP"
2136 call getenv_loc('SCCORPAR',sccorname)
2137 open (isccor,file=sccorname,status='old',action='read')
2138 c print *,"Processor",myrank," opened file ISCCOR"
2139 call getenv_loc('FOURIER',fouriername)
2140 open (ifourier,file=fouriername,status='old',action='read')
2141 c print *,"Processor",myrank," opened file IFOURIER"
2142 call getenv_loc('ELEPAR',elename)
2143 open (ielep,file=elename,status='old',action='read')
2144 c print *,"Processor",myrank," opened file IELEP"
2145 call getenv_loc('SIDEPAR',sidename)
2146 open (isidep,file=sidename,status='old',action='read')
2147 call getenv_loc('LIPTRANPAR',liptranname)
2148 open (iliptranpar,file=liptranname,status='old',action='read')
2149 c print *,"Processor",myrank," opened file ISIDEP"
2150 c print *,"Processor",myrank," opened parameter files"
2152 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2153 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2154 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2155 C Get parameter filenames and open the parameter files.
2156 call getenv_loc('BONDPAR',bondname)
2157 open (ibond,file=bondname,status='old')
2158 call getenv_loc('THETPAR',thetname)
2159 open (ithep,file=thetname,status='old')
2161 call getenv_loc('THETPARPDB',thetname_pdb)
2162 open (ithep_pdb,file=thetname_pdb,status='old')
2164 call getenv_loc('ROTPAR',rotname)
2165 open (irotam,file=rotname,status='old')
2167 call getenv_loc('ROTPARPDB',rotname_pdb)
2168 open (irotam_pdb,file=rotname_pdb,status='old')
2170 call getenv_loc('TORPAR',torname)
2171 open (itorp,file=torname,status='old')
2172 call getenv_loc('TORDPAR',tordname)
2173 open (itordp,file=tordname,status='old')
2174 call getenv_loc('SCCORPAR',sccorname)
2175 open (isccor,file=sccorname,status='old')
2176 call getenv_loc('FOURIER',fouriername)
2177 open (ifourier,file=fouriername,status='old')
2178 call getenv_loc('ELEPAR',elename)
2179 open (ielep,file=elename,status='old')
2180 call getenv_loc('SIDEPAR',sidename)
2181 open (isidep,file=sidename,status='old')
2182 call getenv_loc('LIPTRANPAR',liptranname)
2183 open (iliptranpar,file=liptranname,status='old')
2185 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2187 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2188 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2189 C Get parameter filenames and open the parameter files.
2190 call getenv_loc('BONDPAR',bondname)
2191 open (ibond,file=bondname,status='old',readonly)
2192 call getenv_loc('THETPAR',thetname)
2193 open (ithep,file=thetname,status='old',readonly)
2194 call getenv_loc('ROTPAR',rotname)
2195 open (irotam,file=rotname,status='old',readonly)
2196 call getenv_loc('TORPAR',torname)
2197 open (itorp,file=torname,status='old',readonly)
2198 call getenv_loc('TORDPAR',tordname)
2199 open (itordp,file=tordname,status='old',readonly)
2200 call getenv_loc('SCCORPAR',sccorname)
2201 open (isccor,file=sccorname,status='old',readonly)
2203 call getenv_loc('THETPARPDB',thetname_pdb)
2204 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2206 call getenv_loc('FOURIER',fouriername)
2207 open (ifourier,file=fouriername,status='old',readonly)
2208 call getenv_loc('ELEPAR',elename)
2209 open (ielep,file=elename,status='old',readonly)
2210 call getenv_loc('SIDEPAR',sidename)
2211 open (isidep,file=sidename,status='old',readonly)
2212 call getenv_loc('LIPTRANPAR',liptranname)
2213 open (iliptranpar,file=liptranname,status='old',action='read')
2215 call getenv_loc('ROTPARPDB',rotname_pdb)
2216 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2220 call getenv_loc('TUBEPAR',tubename)
2221 #if defined(WINIFL) || defined(WINPGI)
2222 open (itube,file=tubename,status='old',readonly,shared)
2223 #elif (defined CRAY) || (defined AIX)
2224 open (itube,file=tubename,status='old',action='read')
2226 open (itube,file=tubename,status='old')
2228 open (itube,file=tubename,status='old',readonly)
2233 C 8/9/01 In the newest version SCp interaction constants are read from a file
2234 C Use -DOLDSCP to use hard-coded constants instead.
2236 call getenv_loc('SCPPAR',scpname)
2237 #if defined(WINIFL) || defined(WINPGI)
2238 open (iscpp,file=scpname,status='old',readonly,shared)
2239 #elif (defined CRAY) || (defined AIX)
2240 open (iscpp,file=scpname,status='old',action='read')
2242 open (iscpp,file=scpname,status='old')
2244 open (iscpp,file=scpname,status='old',readonly)
2247 call getenv_loc('PATTERN',patname)
2248 #if defined(WINIFL) || defined(WINPGI)
2249 open (icbase,file=patname,status='old',readonly,shared)
2250 #elif (defined CRAY) || (defined AIX)
2251 open (icbase,file=patname,status='old',action='read')
2253 open (icbase,file=patname,status='old')
2255 open (icbase,file=patname,status='old',readonly)
2258 C Open output file only for CG processes
2259 c print *,"Processor",myrank," fg_rank",fg_rank
2260 if (fg_rank.eq.0) then
2262 if (nodes.eq.1) then
2265 npos = dlog10(dfloat(nodes-1))+1
2267 if (npos.lt.3) npos=3
2268 write (liczba,'(i1)') npos
2269 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2271 write (liczba,form) me
2272 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2273 & liczba(:ilen(liczba))
2274 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2276 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2278 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2279 & liczba(:ilen(liczba))//'.mol2'
2280 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2281 & liczba(:ilen(liczba))//'.stat'
2283 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2284 & //liczba(:ilen(liczba))//'.stat')
2285 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2288 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2289 & liczba(:ilen(liczba))//'.const'
2294 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2295 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2296 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2297 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2298 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2300 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2302 rest2name=prefix(:ilen(prefix))//'.rst'
2304 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2307 #if defined(AIX) || defined(PGI) || defined(CRAY)
2308 if (me.eq.king .or. .not. out1file) then
2309 open(iout,file=outname,status='unknown')
2311 open(iout,file="/dev/null",status="unknown")
2315 if (fg_rank.gt.0) then
2316 write (liczba,'(i3.3)') myrank/nfgtasks
2317 write (ll,'(bz,i3.3)') fg_rank
2318 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2324 open(igeom,file=intname,status='unknown',position='append')
2325 open(ipdb,file=pdbname,status='unknown')
2326 open(imol2,file=mol2name,status='unknown')
2327 open(istat,file=statname,status='unknown',position='append')
2329 c1out open(iout,file=outname,status='unknown')
2332 if (me.eq.king .or. .not.out1file)
2333 & open(iout,file=outname,status='unknown')
2335 if (fg_rank.gt.0) then
2336 write (liczba,'(i3.3)') myrank/nfgtasks
2337 write (ll,'(bz,i3.3)') fg_rank
2338 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2343 open(igeom,file=intname,status='unknown',access='append')
2344 open(ipdb,file=pdbname,status='unknown')
2345 open(imol2,file=mol2name,status='unknown')
2346 open(istat,file=statname,status='unknown',access='append')
2348 c1out open(iout,file=outname,status='unknown')
2351 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2352 csa_seed=prefix(:lenpre)//'.CSA.seed'
2353 csa_history=prefix(:lenpre)//'.CSA.history'
2354 csa_bank=prefix(:lenpre)//'.CSA.bank'
2355 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2356 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2357 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2358 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2359 csa_int=prefix(:lenpre)//'.int'
2360 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2361 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2362 csa_in=prefix(:lenpre)//'.CSA.in'
2363 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2366 write (iout,'(80(1h-))')
2367 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2368 write (iout,'(80(1h-))')
2369 write (iout,*) "Input file : ",
2370 & pref_orig(:ilen(pref_orig))//'.inp'
2371 write (iout,*) "Output file : ",
2372 & outname(:ilen(outname))
2374 write (iout,*) "Sidechain potential file : ",
2375 & sidename(:ilen(sidename))
2377 write (iout,*) "SCp potential file : ",
2378 & scpname(:ilen(scpname))
2380 write (iout,*) "Electrostatic potential file : ",
2381 & elename(:ilen(elename))
2382 write (iout,*) "Cumulant coefficient file : ",
2383 & fouriername(:ilen(fouriername))
2384 write (iout,*) "Torsional parameter file : ",
2385 & torname(:ilen(torname))
2386 write (iout,*) "Double torsional parameter file : ",
2387 & tordname(:ilen(tordname))
2388 write (iout,*) "SCCOR parameter file : ",
2389 & sccorname(:ilen(sccorname))
2390 write (iout,*) "Bond & inertia constant file : ",
2391 & bondname(:ilen(bondname))
2392 write (iout,*) "Bending-torsion parameter file : ",
2393 & thetname(:ilen(thetname))
2394 write (iout,*) "Rotamer parameter file : ",
2395 & rotname(:ilen(rotname))
2396 write (iout,*) "Thetpdb parameter file : ",
2397 & thetname_pdb(:ilen(thetname_pdb))
2398 write (iout,*) "Threading database : ",
2399 & patname(:ilen(patname))
2401 &write (iout,*)" DIRTMP : ",
2403 write (iout,'(80(1h-))')
2407 c----------------------------------------------------------------------------
2408 subroutine card_concat(card)
2409 implicit real*8 (a-h,o-z)
2410 include 'DIMENSIONS'
2411 include 'COMMON.IOUNITS'
2413 character*80 karta,ucase
2415 read (inp,'(a)') karta
2418 do while (karta(80:80).eq.'&')
2419 card=card(:ilen(card)+1)//karta(:79)
2420 read (inp,'(a)') karta
2423 card=card(:ilen(card)+1)//karta
2426 c------------------------------------------------------------------------------
2429 include 'DIMENSIONS'
2430 include 'COMMON.CHAIN'
2431 include 'COMMON.IOUNITS'
2432 include 'COMMON.CONTROL'
2434 include 'COMMON.QRESTR'
2436 include 'COMMON.LAGRANGE.5diag'
2438 include 'COMMON.LAGRANGE'
2441 open(irest2,file=rest2name,status='unknown')
2442 read(irest2,*) totT,EK,potE,totE,t_bath
2445 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2448 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2451 read (irest2,*) iset
2456 c------------------------------------------------------------------------------
2457 subroutine read_fragments
2459 include 'DIMENSIONS'
2463 include 'COMMON.SETUP'
2464 include 'COMMON.CHAIN'
2465 include 'COMMON.IOUNITS'
2467 include 'COMMON.QRESTR'
2468 include 'COMMON.CONTROL'
2470 read(inp,*) nset,nfrag,npair,nfrag_back
2471 loc_qlike=(nfrag_back.lt.0)
2472 nfrag_back=iabs(nfrag_back)
2473 if(me.eq.king.or..not.out1file)
2474 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2475 & " nfrag_back",nfrag_back," loc_qlike",loc_qlike
2477 read(inp,*) mset(iset)
2479 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2481 if(me.eq.king.or..not.out1file)
2482 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2483 & ifrag(2,i,iset), qinfrag(i,iset)
2486 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2488 if(me.eq.king.or..not.out1file)
2489 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2490 & ipair(2,i,iset), qinpair(i,iset)
2494 read(inp,*) wfrag_back(1,i,iset),qin_back(1,i,iset),
2495 & wfrag_back(2,i,iset),qin_back(2,i,iset),
2496 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2497 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2498 if(me.eq.king.or..not.out1file)
2499 & write(iout,*) "A",i,wfrag_back(1,i,iset),qin_back(2,i,iset),
2500 & wfrag_back(2,i,iset),qin_back(3,i,iset),
2501 & wfrag_back(3,i,iset),qin_back(3,i,iset),
2502 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2506 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2507 & wfrag_back(3,i,iset),
2508 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2509 if(me.eq.king.or..not.out1file)
2510 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2511 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2517 C---------------------------------------------------------------------------
2518 subroutine read_afminp
2520 include 'DIMENSIONS'
2524 include 'COMMON.SETUP'
2525 include 'COMMON.CONTROL'
2526 include 'COMMON.CHAIN'
2527 include 'COMMON.IOUNITS'
2528 include 'COMMON.SBRIDGE'
2529 character*320 afmcard
2531 c print *, "wchodze"
2532 call card_concat(afmcard)
2533 call readi(afmcard,"BEG",afmbeg,0)
2534 call readi(afmcard,"END",afmend,0)
2535 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2536 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2537 c print *,'FORCE=' ,forceAFMconst
2538 CCCC NOW PROPERTIES FOR AFM
2541 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2543 distafminit=dsqrt(distafminit)
2544 c print *,'initdist',distafminit
2547 c-------------------------------------------------------------------------------
2548 subroutine read_saxs_constr
2550 include 'DIMENSIONS'
2554 include 'COMMON.SETUP'
2555 include 'COMMON.CONTROL'
2556 include 'COMMON.SAXS'
2557 include 'COMMON.CHAIN'
2558 include 'COMMON.IOUNITS'
2559 include 'COMMON.SBRIDGE'
2560 double precision cm(3),cnorm
2563 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2565 if (saxs_mode.eq.0) then
2566 c SAXS distance distribution
2568 read(inp,*) distsaxs(i),Psaxs(i)
2572 Cnorm = Cnorm + Psaxs(i)
2574 write (iout,*) "Cnorm",Cnorm
2576 Psaxs(i)=Psaxs(i)/Cnorm
2578 write (iout,*) "Normalized distance distribution from SAXS"
2580 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2584 Wsaxs0=Wsaxs0-Psaxs(i)*dlog(Psaxs(i))
2586 write (iout,*) "Wsaxs0",Wsaxs0
2590 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2597 cm(j)=cm(j)+Csaxs(j,i)
2605 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2608 write (iout,*) "SAXS sphere coordinates"
2610 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2616 c-------------------------------------------------------------------------------
2617 subroutine read_dist_constr
2619 include 'DIMENSIONS'
2623 include 'COMMON.SETUP'
2624 include 'COMMON.CONTROL'
2625 include 'COMMON.CHAIN'
2626 include 'COMMON.IOUNITS'
2627 include 'COMMON.SBRIDGE'
2628 include 'COMMON.INTERACT'
2629 integer i,j,k,ii,jj,itemp,link_type,iiend,jjend,kk
2630 integer nfrag_,npair_,ndist_,ifrag_(2,100),ipair_(2,1000)
2631 double precision wfrag_(100),wpair_(1000)
2632 double precision ddjk,dist,dist_cut,fordepthmax
2633 character*5000 controlcard
2634 logical normalize,next
2636 double precision scal_bfac
2637 double precision xlink(4,0:4) /
2639 & 0.0d0,0.0d0,0.0d0,0.0d0, ! default, no xlink potential
2640 & 0.00305218d0,9.46638d0,4.68901d0,4.74347d0, ! ZL
2641 & 0.00214928d0,12.7517d0,0.00375009d0,6.13477d0, ! ADH
2642 & 0.00184547d0,11.2678d0,0.00140292d0,7.00868d0, ! PDH
2643 & 0.000161786d0,6.29273d0,4.40993d0,7.13956d0 / ! DSS
2644 c print *, "WCHODZE"
2645 write (iout,*) "Calling read_dist_constr"
2646 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2648 restr_on_coord=.false.
2657 call card_concat(controlcard)
2658 next = index(controlcard,"NEXT").gt.0
2659 call readi(controlcard,"RESTR_TYPE",restr_type,constr_dist)
2660 write (iout,*) "restr_type",restr_type
2661 call readi(controlcard,"NFRAG",nfrag_,0)
2662 call readi(controlcard,"NFRAG",nfrag_,0)
2663 call readi(controlcard,"NPAIR",npair_,0)
2664 call readi(controlcard,"NDIST",ndist_,0)
2665 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2666 call reada(controlcard,'SCAL_BFAC',scal_bfac,1.0d0)
2667 if (restr_type.eq.10)
2668 & call reada(controlcard,'WBOLTZD',wboltzd,0.591d0)
2669 if (restr_type.eq.12)
2670 & call reada(controlcard,'SCAL_PEAK',scal_peak,5.0d0)
2671 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2672 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2673 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2674 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2675 normalize = index(controlcard,"NORMALIZE").gt.0
2676 write (iout,*) "WBOLTZD",wboltzd
2677 write (iout,*) "SCAL_PEAK",scal_peak
2678 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2679 write (iout,*) "IFRAG"
2681 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2683 write (iout,*) "IPAIR"
2685 write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2687 if (nfrag_.gt.0 .or. restr_type.eq.4 .or. restr_type.eq.5)
2689 & "Distance restraints as generated from reference structure"
2691 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2692 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2693 & ifrag_(2,i)=nstart_sup+nsup-1
2694 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2696 if (wfrag_(i).eq.0.0d0) cycle
2697 do j=ifrag_(1,i),ifrag_(2,i)-1
2698 do k=j+1,ifrag_(2,i)
2699 c write (iout,*) "j",j," k",k
2701 if (restr_type.eq.1) then
2707 forcon(nhpb)=wfrag_(i)
2708 else if (constr_dist.eq.2) then
2709 if (ddjk.le.dist_cut) then
2715 forcon(nhpb)=wfrag_(i)
2717 else if (restr_type.eq.3) then
2723 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2726 if (.not.out1file .or. me.eq.king)
2727 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2728 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2730 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2731 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2737 if (wpair_(i).eq.0.0d0) cycle
2745 do j=ifrag_(1,ii),ifrag_(2,ii)
2746 do k=ifrag_(1,jj),ifrag_(2,jj)
2748 if (restr_type.eq.1) then
2754 forcon(nhpb)=wpair_(i)
2755 else if (constr_dist.eq.2) then
2756 if (ddjk.le.dist_cut) then
2762 forcon(nhpb)=wpair_(i)
2764 else if (restr_type.eq.3) then
2770 forcon(nhpb)=wpair_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2773 if (.not.out1file .or. me.eq.king)
2774 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2775 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2777 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.restr ",
2778 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2785 write (iout,*) "Distance restraints as read from input"
2787 if (restr_type.eq.12) then
2788 read (inp,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2789 & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2790 & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2791 & fordepth_peak(nhpb_peak+1),npeak
2792 c write(iout,*) ihpb_peak(nhpb_peak+1),jhpb_peak(nhpb_peak+1),
2793 c & dhpb_peak(nhpb_peak+1),dhpb1_peak(nhpb_peak+1),
2794 c & ibecarb_peak(nhpb_peak+1),forcon_peak(nhpb_peak+1),
2795 c & fordepth_peak(nhpb_peak+1),npeak
2796 if (forcon_peak(nhpb_peak+1).le.0.0d0.or.
2797 & fordepth_peak(nhpb_peak+1).le.0.0d0)cycle
2798 nhpb_peak=nhpb_peak+1
2799 irestr_type_peak(nhpb_peak)=12
2800 if (ipeak(1,npeak).eq.0) ipeak(1,npeak)=i
2803 if (.not.out1file .or. me.eq.king)
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 write (iout,'(a,5i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2811 & nhpb_peak,ihpb_peak(nhpb_peak),jhpb_peak(nhpb_peak),
2812 & ibecarb_peak(nhpb_peak),npeak,dhpb_peak(nhpb_peak),
2813 & dhpb1_peak(nhpb_peak),forcon_peak(nhpb_peak),
2814 & fordepth_peak(nhpb_peak),irestr_type_peak(nhpb_peak)
2816 if (ibecarb_peak(nhpb_peak).eq.3) then
2817 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2818 else if (ibecarb_peak(nhpb_peak).eq.2) then
2819 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2820 else if (ibecarb_peak(nhpb_peak).eq.1) then
2821 ihpb_peak(nhpb_peak)=ihpb_peak(nhpb_peak)+nres
2822 jhpb_peak(nhpb_peak)=jhpb_peak(nhpb_peak)+nres
2824 else if (restr_type.eq.11) then
2825 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2826 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2827 c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2828 if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle
2830 irestr_type(nhpb)=11
2832 if (.not.out1file .or. me.eq.king)
2833 & write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2834 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2835 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2837 write (iout,'(a,4i5,2f8.2,2f10.5,i5)') "+dist.restr ",
2838 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2839 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),irestr_type(nhpb)
2841 c if (ibecarb(nhpb).gt.0) then
2842 c ihpb(nhpb)=ihpb(nhpb)+nres
2843 c jhpb(nhpb)=jhpb(nhpb)+nres
2845 if (ibecarb(nhpb).eq.3) then
2846 ihpb(nhpb)=ihpb(nhpb)+nres
2847 else if (ibecarb(nhpb).eq.2) then
2848 ihpb(nhpb)=ihpb(nhpb)+nres
2849 else if (ibecarb(nhpb).eq.1) then
2850 ihpb(nhpb)=ihpb(nhpb)+nres
2851 jhpb(nhpb)=jhpb(nhpb)+nres
2853 else if (restr_type.eq.10) then
2854 c Cross-lonk Markov-like potential
2855 call card_concat(controlcard)
2856 call readi(controlcard,"ILINK",ihpb(nhpb+1),0)
2857 call readi(controlcard,"JLINK",jhpb(nhpb+1),0)
2859 if (index(controlcard,"BETA").gt.0) ibecarb(nhpb+1)=1
2860 if (ihpb(nhpb+1).eq.0 .or. jhpb(nhpb+1).eq.0) cycle
2861 if (index(controlcard,"ZL").gt.0) then
2863 else if (index(controlcard,"ADH").gt.0) then
2865 else if (index(controlcard,"PDH").gt.0) then
2867 else if (index(controlcard,"DSS").gt.0) then
2872 call reada(controlcard,"AXLINK",dhpb(nhpb+1),
2873 & xlink(1,link_type))
2874 call reada(controlcard,"BXLINK",dhpb1(nhpb+1),
2875 & xlink(2,link_type))
2876 call reada(controlcard,"CXLINK",fordepth(nhpb+1),
2877 & xlink(3,link_type))
2878 call reada(controlcard,"SIGMA",forcon(nhpb+1),
2879 & xlink(4,link_type))
2880 call reada(controlcard,"SCORE",xlscore(nhpb+1),1.0d0)
2881 c read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),ibecarb(nhpb+1),
2882 c & dhpb(nhpb+1),dhpb1(nhpb+1),forcon(nhpb+1),fordepth(nhpb+1)
2883 if (forcon(nhpb+1).le.0.0d0 .or.
2884 & (dhpb(nhpb+1).eq.0 .and. dhpb1(nhpb+1).eq.0)) cycle
2886 irestr_type(nhpb)=10
2887 if (ibecarb(nhpb).eq.3) then
2888 jhpb(nhpb)=jhpb(nhpb)+nres
2889 else if (ibecarb(nhpb).eq.2) then
2890 ihpb(nhpb)=ihpb(nhpb)+nres
2891 else if (ibecarb(nhpb).eq.1) then
2892 ihpb(nhpb)=ihpb(nhpb)+nres
2893 jhpb(nhpb)=jhpb(nhpb)+nres
2896 if (.not.out1file .or. me.eq.king)
2897 & write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2898 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2899 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2902 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2903 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2904 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2909 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(nhpb+1),
2910 & dhpb1(nhpb+1),ibecarb(nhpb+1),forcon(nhpb+1)
2911 if (forcon(nhpb+1).gt.0.0d0) then
2913 if (dhpb1(nhpb).eq.0.0d0) then
2918 if (ibecarb(nhpb).eq.3) then
2919 jhpb(nhpb)=jhpb(nhpb)+nres
2920 else if (ibecarb(nhpb).eq.2) then
2921 ihpb(nhpb)=ihpb(nhpb)+nres
2922 else if (ibecarb(nhpb).eq.1) then
2923 ihpb(nhpb)=ihpb(nhpb)+nres
2924 jhpb(nhpb)=jhpb(nhpb)+nres
2926 if (dhpb(nhpb).eq.0.0d0)
2927 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2930 if (.not.out1file .or. me.eq.king)
2931 & write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2932 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2934 write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ",
2935 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(nhpb)
2938 C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2939 C if (forcon(nhpb+1).gt.0.0d0) then
2941 C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2944 if (restr_type.eq.4) then
2945 write (iout,*) "The BFAC array"
2947 write (iout,'(i5,f10.5)') i,bfac(i)
2950 if (itype(i).eq.ntyp1) cycle
2952 if (itype(j).eq.ntyp1) cycle
2953 if (itype(i).eq.10) then
2958 if (itype(j).eq.10) then
2968 forcon(nhpb)=scal_bfac**2/(bfac(i)**2+bfac(j)**2)
2971 if (ibecarb(nhpb).gt.0) ibecarb(nhpb)=4-ibecarb(nhpb)
2972 ihpb(nhpb)=i+nres*ii
2973 jhpb(nhpb)=j+nres*jj
2974 dhpb(nhpb)=dist(i+nres*ii,j+nres*jj)
2976 if (.not.out1file .or. me.eq.king) then
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),
2983 write (iout,'(a,4i5,2f8.2,3f10.5,i5)') "+dist.restr ",
2984 & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb),
2985 & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb),xlscore(nhpb),
2995 if (restr_type.eq.5) then
2996 restr_on_coord=.true.
2998 if (itype(i).eq.ntyp1) cycle
2999 bfac(i)=(scal_bfac/bfac(i))**2
3008 if (irestr_type(i).eq.11.and.fordepth(i).gt.fordepthmax)
3009 & fordepthmax=fordepth(i)
3012 if (irestr_type(i).eq.11) fordepth(i)=fordepth(i)/fordepthmax
3017 c-------------------------------------------------------------------------------
3018 subroutine read_constr_homology
3020 include 'DIMENSIONS'
3024 include 'COMMON.SETUP'
3025 include 'COMMON.CONTROL'
3026 include 'COMMON.HOMOLOGY'
3027 include 'COMMON.CHAIN'
3028 include 'COMMON.IOUNITS'
3030 include 'COMMON.QRESTR'
3031 include 'COMMON.GEO'
3032 include 'COMMON.INTERACT'
3033 include 'COMMON.NAMES'
3035 c For new homol impl
3037 include 'COMMON.VAR'
3040 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
3042 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
3043 c & sigma_odl_temp(maxres,maxres,max_template)
3045 character*24 model_ki_dist, model_ki_angle
3046 character*500 controlcard
3047 integer ki,i,ii,j,k,l,ii_in_use(maxdim),i_tmp,idomain_tmp,irec,
3048 & ik,iistart,nres_temp
3051 logical liiflag,lfirst
3054 c FP - Nov. 2014 Temporary specifications for new vars
3056 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
3058 double precision, dimension (max_template,maxres) :: rescore
3059 double precision, dimension (max_template,maxres) :: rescore2
3060 double precision, dimension (max_template,maxres) :: rescore3
3061 double precision distal
3062 character*24 tpl_k_rescore
3063 character*256 pdbfile
3064 c -----------------------------------------------------------------
3065 c Reading multiple PDB ref structures and calculation of retraints
3066 c not using pre-computed ones stored in files model_ki_{dist,angle}
3068 c -----------------------------------------------------------------
3071 c Alternative: reading from input
3072 call card_concat(controlcard)
3073 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
3074 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
3075 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
3076 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
3077 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
3078 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
3079 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
3080 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
3081 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
3082 if(.not.read2sigma.and.start_from_model) then
3083 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
3084 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
3085 start_from_model=.false.
3086 iranconf=(indpdb.le.0)
3088 if(start_from_model .and. (me.eq.king .or. .not. out1file))
3089 & write(iout,*) 'START_FROM_MODELS is ON'
3090 c if(start_from_model .and. rest) then
3091 c if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3092 c write(iout,*) 'START_FROM_MODELS is OFF'
3093 c write(iout,*) 'remove restart keyword from input'
3096 if (start_from_model) nmodel_start=constr_homology
3097 if (homol_nset.gt.1)then
3098 call card_concat(controlcard)
3099 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
3100 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3101 write(iout,*) "iset homology_weight "
3103 write(iout,*) i,waga_homology(i)
3106 iset=mod(kolor,homol_nset)+1
3109 waga_homology(1)=1.0
3112 cd write (iout,*) "nnt",nnt," nct",nct
3119 c write(iout,*) 'nnt=',nnt,'nct=',nct
3122 do k=1,constr_homology
3135 if (read_homol_frag) then
3136 call read_klapaucjusz
3139 do k=1,constr_homology
3141 read(inp,'(a)') pdbfile
3142 pdbfiles_chomo(k)=pdbfile
3143 if(me.eq.king .or. .not. out1file)
3144 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
3145 & pdbfile(:ilen(pdbfile))
3146 open(ipdbin,file=pdbfile,status='old',err=33)
3148 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3149 & pdbfile(:ilen(pdbfile))
3152 c print *,'Begin reading pdb data'
3154 c Files containing res sim or local scores (former containing sigmas)
3157 write(kic2,'(bz,i2.2)') k
3159 tpl_k_rescore="template"//kic2//".sco"
3163 if (read2sigma) then
3164 call readpdb_template(k)
3171 c Distance restraints
3174 C Copy the coordinates from reference coordinates (?)
3175 do i=1,2*nres_chomo(k)
3178 c write (iout,*) "c(",j,i,") =",c(j,i)
3182 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
3184 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
3185 open (ientin,file=tpl_k_rescore,status='old')
3186 if (nnt.gt.1) rescore(k,1)=0.0d0
3187 do irec=nnt,nct ! loop for reading res sim
3188 if (read2sigma) then
3189 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
3190 & rescore3_tmp,idomain_tmp
3192 idomain(k,i_tmp)=idomain_tmp
3193 rescore(k,i_tmp)=rescore_tmp
3194 rescore2(k,i_tmp)=rescore2_tmp
3195 rescore3(k,i_tmp)=rescore3_tmp
3196 if (.not. out1file .or. me.eq.king)
3197 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
3198 & i_tmp,rescore2_tmp,rescore_tmp,
3199 & rescore3_tmp,idomain_tmp
3202 read (ientin,*,end=1401) rescore_tmp
3204 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
3205 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
3206 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
3211 if (waga_dist.ne.0.0d0) then
3219 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3220 c write (iout,*) k,i,j,distal,dist2_cut
3222 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3223 & .and. distal.le.dist2_cut ) then
3229 c write (iout,*) "k",k
3230 c write (iout,*) "i",i," j",j," constr_homology",
3235 if (read2sigma) then
3238 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3240 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3241 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3242 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3244 if (odl(k,ii).le.dist_cut) then
3245 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3248 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3249 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3251 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3252 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3256 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3259 l_homo(k,ii)=.false.
3265 c write (iout,*) "Distance restraints set"
3268 c Theta, dihedral and SC retraints
3270 if (waga_angle.gt.0.0d0) then
3271 c open (ientin,file=tpl_k_sigma_dih,status='old')
3272 c do irec=1,maxres-3 ! loop for reading sigma_dih
3273 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
3274 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
3275 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
3276 c & sigma_dih(k,i+nnt-1)
3281 if (idomain(k,i).eq.0) then
3285 dih(k,i)=phiref(i) ! right?
3286 c read (ientin,*) sigma_dih(k,i) ! original variant
3287 c write (iout,*) "dih(",k,i,") =",dih(k,i)
3288 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3289 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3290 c & "rescore(",k,i-2,") =",rescore(k,i-2),
3291 c & "rescore(",k,i-3,") =",rescore(k,i-3)
3293 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3294 & rescore(k,i-2)+rescore(k,i-3))/4.0
3295 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
3296 c write (iout,*) "Raw sigmas for dihedral angle restraints"
3297 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
3298 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3299 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
3300 c Instead of res sim other local measure of b/b str reliability possible
3301 if (sigma_dih(k,i).ne.0)
3302 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3303 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
3307 c write (iout,*) "Dihedral angle restraints set"
3310 if (waga_theta.gt.0.0d0) then
3311 c open (ientin,file=tpl_k_sigma_theta,status='old')
3312 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
3313 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
3314 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
3315 c & sigma_theta(k,i+nnt-1)
3320 do i = nnt+2,nct ! right? without parallel.
3321 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3322 c do i=ithet_start,ithet_end ! with FG parallel.
3323 if (idomain(k,i).eq.0) then
3324 sigma_theta(k,i)=0.0
3327 thetatpl(k,i)=thetaref(i)
3328 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3329 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3330 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3331 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3332 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3333 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3334 & rescore(k,i-2))/3.0
3335 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3336 if (sigma_theta(k,i).ne.0)
3337 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3339 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3340 c rescore(k,i-2) ! right expression ?
3341 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3344 c write (iout,*) "Angle restraints set"
3347 if (waga_d.gt.0.0d0) then
3348 c open (ientin,file=tpl_k_sigma_d,status='old')
3349 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3350 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3351 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3352 c & sigma_d(k,i+nnt-1)
3356 do i = nnt,nct ! right? without parallel.
3357 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3358 c do i=loc_start,loc_end ! with FG parallel.
3359 if (itype(i).eq.10) cycle
3360 if (idomain(k,i).eq.0 ) then
3367 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3368 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3369 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3370 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3371 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3372 if (sigma_d(k,i).ne.0)
3373 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3375 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3376 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3377 c read (ientin,*) sigma_d(k,i) ! 1st variant
3381 c write (iout,*) "SC restraints set"
3384 c remove distance restraints not used in any model from the list
3385 c shift data in all arrays
3387 c write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
3388 if (waga_dist.ne.0.0d0) then
3395 c if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3396 c & .and. distal.le.dist2_cut ) then
3397 c write (iout,*) "i",i," j",j," ii",ii
3399 if (ii_in_use(ii).eq.0.and.liiflag.or.
3400 & ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
3407 if(i10.eq.lim_odl) i10=i10+1
3409 ires_homo(iistart+ki)=ires_homo(ki+i01)
3410 jres_homo(iistart+ki)=jres_homo(ki+i01)
3411 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
3412 do k=1,constr_homology
3413 odl(k,iistart+ki)=odl(k,ki+i01)
3414 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
3415 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
3418 iistart=iistart+i10-i01
3421 if (ii_in_use(ii).ne.0.and..not.liiflag) then
3429 c write (iout,*) "Removing distances completed"
3431 endif ! .not. klapaucjusz
3433 if (constr_homology.gt.0) call homology_partition
3434 c write (iout,*) "After homology_partition"
3436 if (constr_homology.gt.0) call init_int_table
3437 c write (iout,*) "After init_int_table"
3439 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3440 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3444 if (.not.out_template_restr) return
3445 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3446 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3447 write (iout,*) "Distance restraints from templates"
3449 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3450 & ii,ires_homo(ii),jres_homo(ii),
3451 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3452 & ki=1,constr_homology)
3454 write (iout,*) "Dihedral angle restraints from templates"
3456 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3457 & (rad2deg*dih(ki,i),
3458 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3460 write (iout,*) "Virtual-bond angle restraints from templates"
3462 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3463 & (rad2deg*thetatpl(ki,i),
3464 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3466 write (iout,*) "SC restraints from templates"
3468 write(iout,'(i5,100(4f8.2,4x))') i,
3469 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3470 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3473 c -----------------------------------------------------------------
3476 c----------------------------------------------------------------------
3478 subroutine flush(iu)
3483 subroutine flush(iu)
3488 c------------------------------------------------------------------------------
3489 subroutine copy_to_tmp(source)
3491 include "DIMENSIONS"
3492 include "COMMON.IOUNITS"
3493 character*(*) source
3494 character* 256 tmpfile
3498 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3499 inquire(file=tmpfile,exist=ex)
3501 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3502 & " to temporary directory..."
3503 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3504 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3508 c------------------------------------------------------------------------------
3509 subroutine move_from_tmp(source)
3511 include "DIMENSIONS"
3512 include "COMMON.IOUNITS"
3513 character*(*) source
3516 write (*,*) "Moving ",source(:ilen(source)),
3517 & " from temporary directory to working directory"
3518 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3519 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3522 c------------------------------------------------------------------------------
3523 subroutine random_init(seed)
3525 C Initialize random number generator
3528 include 'DIMENSIONS'
3531 logical OKRandom, prng_restart
3533 integer iseed_array(4)
3534 integer error_msg,ierr
3536 include 'COMMON.IOUNITS'
3537 include 'COMMON.TIME1'
3538 include 'COMMON.THREAD'
3539 include 'COMMON.SBRIDGE'
3540 include 'COMMON.CONTROL'
3541 include 'COMMON.MCM'
3542 include 'COMMON.MAP'
3543 include 'COMMON.HEADER'
3544 include 'COMMON.CSA'
3545 include 'COMMON.CHAIN'
3546 include 'COMMON.MUCA'
3548 include 'COMMON.FFIELD'
3549 include 'COMMON.SETUP'
3551 double precision seed,ran_number
3552 iseed=-dint(dabs(seed))
3553 if (iseed.eq.0) then
3554 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3555 & 'Random seed undefined. The program will stop.'
3556 write (*,'(/80(1h*)/20x,a/80(1h*))')
3557 & 'Random seed undefined. The program will stop.'
3559 call mpi_finalize(mpi_comm_world,ierr)
3561 stop 'Bad random seed.'
3564 if (fg_rank.eq.0) then
3568 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
3569 OKRandom = prng_restart(me,iseed)
3572 tmp=65536.0d0**(4-i)
3573 iseed_array(i) = dint(seed/tmp)
3574 seed=seed-iseed_array(i)*tmp
3577 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3578 & (iseed_array(i),i=1,4)
3579 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3580 & (iseed_array(i),i=1,4)
3581 OKRandom = prng_restart(me,iseed_array)
3584 c r1 = prng_next(me)
3585 r1=ran_number(0.0D0,1.0D0)
3587 & write (iout,*) 'ran_num',r1
3588 if (r1.lt.0.0d0) OKRandom=.false.
3590 if (.not.OKRandom) then
3591 write (iout,*) 'PRNG IS NOT WORKING!!!'
3592 print *,'PRNG IS NOT WORKING!!!'
3595 call mpi_abort(mpi_comm_world,error_msg,ierr)
3598 write (iout,*) 'too many processors for parallel prng'
3599 write (*,*) 'too many processors for parallel prng'
3607 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3611 c----------------------------------------------------------------------
3612 subroutine read_klapaucjusz
3614 include 'DIMENSIONS'
3618 include 'COMMON.SETUP'
3619 include 'COMMON.CONTROL'
3620 include 'COMMON.HOMOLOGY'
3621 include 'COMMON.CHAIN'
3622 include 'COMMON.IOUNITS'
3624 include 'COMMON.GEO'
3625 include 'COMMON.INTERACT'
3626 include 'COMMON.NAMES'
3627 character*256 fragfile
3628 integer ninclust(maxclust),inclust(max_template,maxclust),
3629 & nresclust(maxclust),iresclust(maxres,maxclust),nclust
3632 character*24 model_ki_dist, model_ki_angle
3633 character*500 controlcard
3634 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp,
3635 & ik,ll,ii,kk,iistart,iishift,lim_xx
3636 double precision distal
3637 logical lprn /.true./
3644 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
3645 double precision, dimension (max_template,maxres) :: rescore
3646 double precision, dimension (max_template,maxres) :: rescore2
3647 character*24 tpl_k_rescore
3648 character*256 pdbfile
3651 c For new homol impl
3653 include 'COMMON.VAR'
3655 call getenv("FRAGFILE",fragfile)
3656 open(ientin,file=fragfile,status="old",err=10)
3657 read(ientin,*) constr_homology,nclust
3658 nmodel_start=constr_homology
3664 do k=1,constr_homology
3665 read(ientin,'(a)') pdbfile
3666 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file',
3667 & pdbfile(:ilen(pdbfile))
3668 pdbfiles_chomo(k)=pdbfile
3669 open(ipdbin,file=pdbfile,status='old',err=33)
3671 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
3672 & pdbfile(:ilen(pdbfile))
3677 call readpdb_template(k)
3687 read(ientin,*) ninclust(i),nresclust(i)
3688 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
3689 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
3692 c Loop over clusters
3695 do ll = 1,ninclust(l)
3703 idomain(k,iresclust(i,l)+1) = 1
3705 idomain(k,iresclust(i,l)) = 1
3709 c Distance restraints
3712 C Copy the coordinates from reference coordinates (?)
3718 c write (iout,*) "c(",j,i,") =",c(j,i)
3721 call int_from_cart(.true.,.false.)
3722 call sc_loc_geom(.false.)
3724 thetaref(i)=theta(i)
3728 if (waga_dist.ne.0.0d0) then
3736 distal=dsqrt(x12*x12+y12*y12+z12*z12)
3737 c write (iout,*) k,i,j,distal,dist2_cut
3739 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
3740 & .and. distal.le.dist2_cut ) then
3746 c write (iout,*) "k",k
3747 c write (iout,*) "i",i," j",j," constr_homology",
3752 if (read2sigma) then
3755 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
3757 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
3758 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
3759 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3761 if (odl(k,ii).le.dist_cut) then
3762 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
3765 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3766 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
3768 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
3769 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
3773 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
3776 c l_homo(k,ii)=.false.
3783 c Theta, dihedral and SC retraints
3785 if (waga_angle.gt.0.0d0) then
3787 if (idomain(k,i).eq.0) then
3788 c sigma_dih(k,i)=0.0
3792 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
3793 & rescore(k,i-2)+rescore(k,i-3))/4.0
3794 c write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
3795 c & " sigma_dihed",sigma_dih(k,i)
3796 if (sigma_dih(k,i).ne.0)
3797 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
3802 if (waga_theta.gt.0.0d0) then
3804 if (idomain(k,i).eq.0) then
3805 c sigma_theta(k,i)=0.0
3808 thetatpl(k,i)=thetaref(i)
3809 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3810 & rescore(k,i-2))/3.0
3811 if (sigma_theta(k,i).ne.0)
3812 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3816 if (waga_d.gt.0.0d0) then
3818 if (itype(i).eq.10) cycle
3819 if (idomain(k,i).eq.0 ) then
3826 sigma_d(k,i)=rescore(k,i)
3827 if (sigma_d(k,i).ne.0)
3828 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3829 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
3835 c remove distance restraints not used in any model from the list
3836 c shift data in all arrays
3838 if (waga_dist.ne.0.0d0) then
3844 if (ii_in_use(ii).eq.0.and.liiflag) then
3848 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3849 & .not.liiflag.and.ii.eq.lim_odl) then
3850 if (ii.eq.lim_odl) then
3851 iishift=ii-iistart+1
3856 do ki=iistart,lim_odl-iishift
3857 ires_homo(ki)=ires_homo(ki+iishift)
3858 jres_homo(ki)=jres_homo(ki+iishift)
3859 ii_in_use(ki)=ii_in_use(ki+iishift)
3860 do k=1,constr_homology
3861 odl(k,ki)=odl(k,ki+iishift)
3862 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3863 l_homo(k,ki)=l_homo(k,ki+iishift)
3867 lim_odl=lim_odl-iishift
3874 10 stop "Error in fragment file"