8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
743 integer :: maxinter,junk,kk,ii,ncatprotparm
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2
750 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
751 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
752 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 ! For printing parameters after they are read set the following in the UNRES
757 ! setenv PRINT_PARM YES
759 ! To print parameters in LaTeX format rather than as ASCII tables:
763 call getenv_loc("PRINT_PARM",lancuch)
764 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
765 call getenv_loc("LATEX",lancuch)
766 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
768 dwa16=2.0d0**(1.0d0/6.0d0)
770 ! Assign virtual-bond length
773 vblinv2=vblinv*vblinv
775 ! Read the virtual-bond parameters, masses, and moments of inertia
776 ! and Stokes' radii of the peptide group and side chains
778 allocate(dsc(ntyp1)) !(ntyp1)
779 allocate(dsc_inv(ntyp1)) !(ntyp1)
780 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
781 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
782 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
783 allocate(nbondterm(ntyp)) !(ntyp)
784 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
785 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
786 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(long_r_sidechain(ntyp))
788 allocate(short_r_sidechain(ntyp))
793 allocate(msc(ntyp+1)) !(ntyp+1)
794 allocate(isc(ntyp+1)) !(ntyp+1)
795 allocate(restok(ntyp+1)) !(ntyp+1)
797 read (ibond,*) vbldp0,akp,mp,ip,pstok
800 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
801 dsc(i) = vbldsc0(1,i)
805 dsc_inv(i)=1.0D0/dsc(i)
814 allocate(msc(ntyp+1,5)) !(ntyp+1)
815 allocate(isc(ntyp+1,5)) !(ntyp+1)
816 allocate(restok(ntyp+1,5)) !(ntyp+1)
818 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
820 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
821 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
822 dsc(i) = vbldsc0(1,i)
826 dsc_inv(i)=1.0D0/dsc(i)
830 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
833 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
834 ! dsc(i) = vbldsc0_nucl(1,i)
838 ! dsc_inv(i)=1.0D0/dsc(i)
841 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
842 ! do i=1,ntyp_molec(2)
843 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
844 ! aksc_nucl(j,i),abond0_nucl(j,i),&
845 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
846 ! dsc(i) = vbldsc0(1,i)
850 ! dsc_inv(i)=1.0D0/dsc(i)
855 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
856 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
858 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
860 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
861 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
863 write (iout,'(13x,3f10.5)') &
864 vbldsc0(j,i),aksc(j,i),abond0(j,i)
869 read(iion,*) msc(i,5),restok(i,5)
870 print *,msc(i,5),restok(i,5)
874 read (iion,*) ncatprotparm
875 allocate(catprm(ncatprotparm,4))
877 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
880 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
881 !----------------------------------------------------
882 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
883 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
884 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
885 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
886 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
887 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
897 allocate(liptranene(ntyp))
898 !C reading lipid parameters
899 write (iout,*) "iliptranpar",iliptranpar
901 read(iliptranpar,*) pepliptran
904 read(iliptranpar,*) liptranene(i)
905 print *,liptranene(i)
911 ! Read the parameters of the probability distribution/energy expression
912 ! of the virtual-bond valence angles theta
915 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
916 (bthet(j,i,1,1),j=1,2)
917 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
918 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
919 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
923 athet(1,i,1,-1)=athet(1,i,1,1)
924 athet(2,i,1,-1)=athet(2,i,1,1)
925 bthet(1,i,1,-1)=-bthet(1,i,1,1)
926 bthet(2,i,1,-1)=-bthet(2,i,1,1)
927 athet(1,i,-1,1)=-athet(1,i,1,1)
928 athet(2,i,-1,1)=-athet(2,i,1,1)
929 bthet(1,i,-1,1)=bthet(1,i,1,1)
930 bthet(2,i,-1,1)=bthet(2,i,1,1)
934 athet(1,i,-1,-1)=athet(1,-i,1,1)
935 athet(2,i,-1,-1)=-athet(2,-i,1,1)
936 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
937 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
938 athet(1,i,-1,1)=athet(1,-i,1,1)
939 athet(2,i,-1,1)=-athet(2,-i,1,1)
940 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
941 bthet(2,i,-1,1)=bthet(2,-i,1,1)
942 athet(1,i,1,-1)=-athet(1,-i,1,1)
943 athet(2,i,1,-1)=athet(2,-i,1,1)
944 bthet(1,i,1,-1)=bthet(1,-i,1,1)
945 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
950 polthet(j,i)=polthet(j,-i)
953 gthet(j,i)=gthet(j,-i)
961 'Parameters of the virtual-bond valence angles:'
962 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
963 ' ATHETA0 ',' A1 ',' A2 ',&
966 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
967 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
969 write (iout,'(/a/9x,5a/79(1h-))') &
970 'Parameters of the expression for sigma(theta_c):',&
971 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
972 ' ALPH3 ',' SIGMA0C '
974 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
975 (polthet(j,i),j=0,3),sigc0(i)
977 write (iout,'(/a/9x,5a/79(1h-))') &
978 'Parameters of the second gaussian:',&
979 ' THETA0 ',' SIGMA0 ',' G1 ',&
982 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
983 sig0(i),(gthet(j,i),j=1,3)
987 'Parameters of the virtual-bond valence angles:'
988 write (iout,'(/a/9x,5a/79(1h-))') &
989 'Coefficients of expansion',&
990 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
991 ' b1*10^1 ',' b2*10^1 '
993 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
994 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
995 (10*bthet(j,i,1,1),j=1,2)
997 write (iout,'(/a/9x,5a/79(1h-))') &
998 'Parameters of the expression for sigma(theta_c):',&
999 ' alpha0 ',' alph1 ',' alph2 ',&
1000 ' alhp3 ',' sigma0c '
1002 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1003 (polthet(j,i),j=0,3),sigc0(i)
1005 write (iout,'(/a/9x,5a/79(1h-))') &
1006 'Parameters of the second gaussian:',&
1007 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1010 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1011 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1017 ! Read the parameters of Utheta determined from ab initio surfaces
1018 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1020 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1021 ntheterm3,nsingle,ndouble
1022 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1024 !----------------------------------------------------
1025 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1026 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1027 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1028 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1029 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1030 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1031 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1032 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1033 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1034 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1035 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1036 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1037 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1038 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1039 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1040 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1041 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1042 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1043 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1044 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1045 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1046 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1047 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1048 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1050 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1052 ithetyp(i)=-ithetyp(-i)
1055 aa0thet(:,:,:,:)=0.0d0
1056 aathet(:,:,:,:,:)=0.0d0
1057 bbthet(:,:,:,:,:,:)=0.0d0
1058 ccthet(:,:,:,:,:,:)=0.0d0
1059 ddthet(:,:,:,:,:,:)=0.0d0
1060 eethet(:,:,:,:,:,:)=0.0d0
1061 ffthet(:,:,:,:,:,:,:)=0.0d0
1062 ggthet(:,:,:,:,:,:,:)=0.0d0
1064 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1066 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1067 ! VAR:1=non-glicyne non-proline 2=proline
1068 ! VAR:negative values for D-aminoacid
1070 do j=-nthetyp,nthetyp
1071 do k=-nthetyp,nthetyp
1072 read (ithep,'(6a)',end=111,err=111) res1
1073 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1074 ! VAR: aa0thet is variable describing the average value of Foureir
1075 ! VAR: expansion series
1076 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1077 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1078 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1079 read (ithep,*,end=111,err=111) &
1080 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1081 read (ithep,*,end=111,err=111) &
1082 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1083 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1084 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1085 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1087 read (ithep,*,end=111,err=111) &
1088 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1089 ffthet(lll,llll,ll,i,j,k,iblock),&
1090 ggthet(llll,lll,ll,i,j,k,iblock),&
1091 ggthet(lll,llll,ll,i,j,k,iblock),&
1092 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1097 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1098 ! coefficients of theta-and-gamma-dependent terms are zero.
1099 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1100 ! RECOMENTDED AFTER VERSION 3.3)
1104 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1105 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1107 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1108 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1111 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1113 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1116 ! AND COMMENT THE LOOPS BELOW
1120 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1121 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1123 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1124 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1127 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1129 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1134 ! Substitution for D aminoacids from symmetry.
1137 do j=-nthetyp,nthetyp
1138 do k=-nthetyp,nthetyp
1139 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1141 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1145 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1146 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1147 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1148 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1154 ffthet(llll,lll,ll,i,j,k,iblock)= &
1155 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1156 ffthet(lll,llll,ll,i,j,k,iblock)= &
1157 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1158 ggthet(llll,lll,ll,i,j,k,iblock)= &
1159 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1160 ggthet(lll,llll,ll,i,j,k,iblock)= &
1161 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1170 ! Control printout of the coefficients of virtual-bond-angle potentials
1173 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1178 write (iout,'(//4a)') &
1179 'Type ',onelett(i),onelett(j),onelett(k)
1180 write (iout,'(//a,10x,a)') " l","a[l]"
1181 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1182 write (iout,'(i2,1pe15.5)') &
1183 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1185 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1186 "b",l,"c",l,"d",l,"e",l
1188 write (iout,'(i2,4(1pe15.5))') m,&
1189 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1190 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1194 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1195 "f+",l,"f-",l,"g+",l,"g-",l
1198 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1199 ffthet(n,m,l,i,j,k,iblock),&
1200 ffthet(m,n,l,i,j,k,iblock),&
1201 ggthet(n,m,l,i,j,k,iblock),&
1202 ggthet(m,n,l,i,j,k,iblock)
1212 write (2,*) "Start reading THETA_PDB",ithep_pdb
1214 ! write (2,*) 'i=',i
1215 read (ithep_pdb,*,err=111,end=111) &
1216 a0thet(i),(athet(j,i,1,1),j=1,2),&
1217 (bthet(j,i,1,1),j=1,2)
1218 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1219 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1220 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1221 sigc0(i)=sigc0(i)**2
1224 athet(1,i,1,-1)=athet(1,i,1,1)
1225 athet(2,i,1,-1)=athet(2,i,1,1)
1226 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1227 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1228 athet(1,i,-1,1)=-athet(1,i,1,1)
1229 athet(2,i,-1,1)=-athet(2,i,1,1)
1230 bthet(1,i,-1,1)=bthet(1,i,1,1)
1231 bthet(2,i,-1,1)=bthet(2,i,1,1)
1234 a0thet(i)=a0thet(-i)
1235 athet(1,i,-1,-1)=athet(1,-i,1,1)
1236 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1237 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1238 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1239 athet(1,i,-1,1)=athet(1,-i,1,1)
1240 athet(2,i,-1,1)=-athet(2,-i,1,1)
1241 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1242 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1243 athet(1,i,1,-1)=-athet(1,-i,1,1)
1244 athet(2,i,1,-1)=athet(2,-i,1,1)
1245 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1246 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1247 theta0(i)=theta0(-i)
1251 polthet(j,i)=polthet(j,-i)
1254 gthet(j,i)=gthet(j,-i)
1257 write (2,*) "End reading THETA_PDB"
1261 !--------------- Reading theta parameters for nucleic acid-------
1262 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1263 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1264 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1265 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1266 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1267 nthetyp_nucl+1,nthetyp_nucl+1))
1268 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1269 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1270 nthetyp_nucl+1,nthetyp_nucl+1))
1271 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1272 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1273 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1274 nthetyp_nucl+1,nthetyp_nucl+1))
1275 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1276 nthetyp_nucl+1,nthetyp_nucl+1))
1277 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1278 nthetyp_nucl+1,nthetyp_nucl+1))
1279 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1280 nthetyp_nucl+1,nthetyp_nucl+1))
1281 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1282 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1283 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1284 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1285 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1286 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1288 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1289 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1291 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1293 aa0thet_nucl(:,:,:)=0.0d0
1294 aathet_nucl(:,:,:,:)=0.0d0
1295 bbthet_nucl(:,:,:,:,:)=0.0d0
1296 ccthet_nucl(:,:,:,:,:)=0.0d0
1297 ddthet_nucl(:,:,:,:,:)=0.0d0
1298 eethet_nucl(:,:,:,:,:)=0.0d0
1299 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1300 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1305 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1306 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1307 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1308 read (ithep_nucl,*,end=111,err=111) &
1309 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1310 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1311 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1312 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1313 read (ithep_nucl,*,end=111,err=111) &
1314 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1315 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1316 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1321 !-------------------------------------------
1322 allocate(nlob(ntyp1)) !(ntyp1)
1323 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1324 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1325 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1332 gaussc(:,:,:,:)=0.0D0
1336 ! Read the parameters of the probability distribution/energy expression
1337 ! of the side chains.
1340 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1344 dsc_inv(i)=1.0D0/dsc(i)
1355 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1356 ((blower(k,l,1),l=1,k),k=1,3)
1357 censc(1,1,-i)=censc(1,1,i)
1358 censc(2,1,-i)=censc(2,1,i)
1359 censc(3,1,-i)=-censc(3,1,i)
1361 read (irotam,*,end=112,err=112) bsc(j,i)
1362 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1363 ((blower(k,l,j),l=1,k),k=1,3)
1364 censc(1,j,-i)=censc(1,j,i)
1365 censc(2,j,-i)=censc(2,j,i)
1366 censc(3,j,-i)=-censc(3,j,i)
1367 ! BSC is amplitude of Gaussian
1374 akl=akl+blower(k,m,j)*blower(l,m,j)
1378 if (((k.eq.3).and.(l.ne.3)) &
1379 .or.((l.eq.3).and.(k.ne.3))) then
1380 gaussc(k,l,j,-i)=-akl
1381 gaussc(l,k,j,-i)=-akl
1383 gaussc(k,l,j,-i)=akl
1384 gaussc(l,k,j,-i)=akl
1393 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1396 if (nlobi.gt.0) then
1398 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1399 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1400 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1401 'log h',(bsc(j,i),j=1,nlobi)
1402 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1403 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1405 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1406 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1409 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1410 write (iout,'(a,f10.4,4(16x,f10.4))') &
1411 'Center ',(bsc(j,i),j=1,nlobi)
1412 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1421 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1422 ! added by Urszula Kozlowska 07/11/2007
1424 !el Maximum number of SC local term fitting function coefficiants
1425 !el integer,parameter :: maxsccoef=65
1427 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1430 read (irotam,*,end=112,err=112)
1432 read (irotam,*,end=112,err=112)
1435 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1439 !---------reading nucleic acid parameters for rotamers-------------------
1440 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1441 do i=1,ntyp_molec(2)
1442 read (irotam_nucl,*,end=112,err=112)
1444 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1450 write (iout,*) "Base rotamer parameters"
1451 do i=1,ntyp_molec(2)
1452 write (iout,'(a)') restyp(i,2)
1453 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1458 ! Read the parameters of the probability distribution/energy expression
1459 ! of the side chains.
1461 write (2,*) "Start reading ROTAM_PDB"
1463 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1467 dsc_inv(i)=1.0D0/dsc(i)
1478 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1479 ((blower(k,l,1),l=1,k),k=1,3)
1481 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1482 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1483 ((blower(k,l,j),l=1,k),k=1,3)
1490 akl=akl+blower(k,m,j)*blower(l,m,j)
1500 write (2,*) "End reading ROTAM_PDB"
1506 ! Read torsional parameters in old format
1508 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1510 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1511 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1512 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1514 !el from energy module--------
1515 allocate(v1(nterm_old,ntortyp,ntortyp))
1516 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1517 !el---------------------------
1522 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1528 write (iout,'(/a/)') 'Torsional constants:'
1531 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1532 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1538 ! Read torsional parameters
1540 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1541 read (itorp,*,end=113,err=113) ntortyp
1542 !el from energy module---------
1543 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1544 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1546 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1547 allocate(vlor2(maxlor,ntortyp,ntortyp))
1548 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1549 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1551 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1552 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1553 !el---------------------------
1556 !el---------------------------
1558 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1560 itortyp(i)=-itortyp(-i)
1562 itortyp(ntyp1)=ntortyp
1563 itortyp(-ntyp1)=-ntortyp
1565 write (iout,*) 'ntortyp',ntortyp
1567 do j=-ntortyp+1,ntortyp-1
1568 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1570 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1571 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1574 do k=1,nterm(i,j,iblock)
1575 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1577 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1578 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1579 v0ij=v0ij+si*v1(k,i,j,iblock)
1581 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1582 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1583 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1585 do k=1,nlor(i,j,iblock)
1586 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1587 vlor2(k,i,j),vlor3(k,i,j)
1588 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1591 v0(-i,-j,iblock)=v0ij
1597 write (iout,'(/a/)') 'Torsional constants:'
1599 do i=-ntortyp,ntortyp
1600 do j=-ntortyp,ntortyp
1601 write (iout,*) 'ityp',i,' jtyp',j
1602 write (iout,*) 'Fourier constants'
1603 do k=1,nterm(i,j,iblock)
1604 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1607 write (iout,*) 'Lorenz constants'
1608 do k=1,nlor(i,j,iblock)
1609 write (iout,'(3(1pe15.5))') &
1610 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1616 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1618 ! 6/23/01 Read parameters for double torsionals
1620 !el from energy module------------
1621 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1622 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1623 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1624 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1625 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1626 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1627 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1628 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1629 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1630 !---------------------------------
1634 do j=-ntortyp+1,ntortyp-1
1635 do k=-ntortyp+1,ntortyp-1
1636 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1637 ! write (iout,*) "OK onelett",
1640 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1641 .or. t3.ne.toronelet(k)) then
1642 write (iout,*) "Error in double torsional parameter file",&
1645 call MPI_Finalize(Ierror)
1647 stop "Error in double torsional parameter file"
1649 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1650 ntermd_2(i,j,k,iblock)
1651 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1652 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1653 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1654 ntermd_1(i,j,k,iblock))
1655 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1656 ntermd_1(i,j,k,iblock))
1657 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1658 ntermd_1(i,j,k,iblock))
1659 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1660 ntermd_1(i,j,k,iblock))
1661 ! Martix of D parameters for one dimesional foureir series
1662 do l=1,ntermd_1(i,j,k,iblock)
1663 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1664 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1665 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1666 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1667 ! write(iout,*) "whcodze" ,
1668 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1670 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1671 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1672 v2s(m,l,i,j,k,iblock),&
1673 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1674 ! Martix of D parameters for two dimesional fourier series
1675 do l=1,ntermd_2(i,j,k,iblock)
1677 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1678 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1679 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1680 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1689 write (iout,*) 'Constants for double torsionals'
1692 do j=-ntortyp+1,ntortyp-1
1693 do k=-ntortyp+1,ntortyp-1
1694 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1695 ' nsingle',ntermd_1(i,j,k,iblock),&
1696 ' ndouble',ntermd_2(i,j,k,iblock)
1698 write (iout,*) 'Single angles:'
1699 do l=1,ntermd_1(i,j,k,iblock)
1700 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1701 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1702 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1703 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1706 write (iout,*) 'Pairs of angles:'
1707 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1708 do l=1,ntermd_2(i,j,k,iblock)
1709 write (iout,'(i5,20f10.5)') &
1710 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1713 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1714 do l=1,ntermd_2(i,j,k,iblock)
1715 write (iout,'(i5,20f10.5)') &
1716 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1717 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1726 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1727 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
1728 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
1729 !el from energy module---------
1730 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1731 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1733 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1734 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
1735 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
1736 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1738 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
1739 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1740 !el---------------------------
1743 !el--------------------
1744 read (itorp_nucl,*,end=113,err=113) &
1745 (itortyp_nucl(i),i=1,ntyp_molec(2))
1746 ! print *,itortyp_nucl(:)
1747 !c write (iout,*) 'ntortyp',ntortyp
1750 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
1751 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
1754 do k=1,nterm_nucl(i,j)
1755 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
1756 v0ij=v0ij+si*v1_nucl(k,i,j)
1759 do k=1,nlor_nucl(i,j)
1760 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
1761 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
1762 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
1768 ! Read of Side-chain backbone correlation parameters
1769 ! Modified 11 May 2012 by Adasko
1772 read (isccor,*,end=119,err=119) nsccortyp
1774 !el from module energy-------------
1775 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1776 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1777 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1778 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1779 !-----------------------------------
1781 !el from module energy-------------
1782 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1784 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1786 isccortyp(i)=-isccortyp(-i)
1788 iscprol=isccortyp(20)
1789 ! write (iout,*) 'ntortyp',ntortyp
1791 !c maxinter is maximum interaction sites
1792 !el from module energy---------
1793 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1794 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1795 -nsccortyp:nsccortyp))
1796 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1797 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1798 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1799 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1800 !-----------------------------------
1804 read (isccor,*,end=119,err=119) &
1805 nterm_sccor(i,j),nlor_sccor(i,j)
1811 nterm_sccor(-i,j)=nterm_sccor(i,j)
1812 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1813 nterm_sccor(i,-j)=nterm_sccor(i,j)
1814 do k=1,nterm_sccor(i,j)
1815 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1817 if (j.eq.iscprol) then
1818 if (i.eq.isccortyp(10)) then
1819 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1820 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1822 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1823 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1824 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1825 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1826 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1827 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1828 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1829 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1832 if (i.eq.isccortyp(10)) then
1833 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1834 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1836 if (j.eq.isccortyp(10)) then
1837 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1838 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1840 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1841 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1842 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1843 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1844 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1845 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1849 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1850 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1851 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1852 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1855 do k=1,nlor_sccor(i,j)
1856 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1857 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1858 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1859 (1+vlor3sccor(k,i,j)**2)
1861 v0sccor(l,i,j)=v0ijsccor
1862 v0sccor(l,-i,j)=v0ijsccor1
1863 v0sccor(l,i,-j)=v0ijsccor2
1864 v0sccor(l,-i,-j)=v0ijsccor3
1870 !el from module energy-------------
1871 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1873 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1874 ! write (iout,*) 'ntortyp',ntortyp
1876 !c maxinter is maximum interaction sites
1877 !el from module energy---------
1878 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1879 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1880 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1881 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1882 !-----------------------------------
1886 read (isccor,*,end=119,err=119) &
1887 nterm_sccor(i,j),nlor_sccor(i,j)
1891 do k=1,nterm_sccor(i,j)
1892 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1894 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1897 do k=1,nlor_sccor(i,j)
1898 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1899 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1900 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1901 (1+vlor3sccor(k,i,j)**2)
1903 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1912 write (iout,'(/a/)') 'Torsional constants:'
1915 write (iout,*) 'ityp',i,' jtyp',j
1916 write (iout,*) 'Fourier constants'
1917 do k=1,nterm_sccor(i,j)
1918 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1920 write (iout,*) 'Lorenz constants'
1921 do k=1,nlor_sccor(i,j)
1922 write (iout,'(3(1pe15.5))') &
1923 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1930 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1931 ! interaction energy of the Gly, Ala, and Pro prototypes.
1936 write (iout,*) "Coefficients of the cumulants"
1938 read (ifourier,*) nloctyp
1939 !write(iout,*) "nloctyp",nloctyp
1940 !el from module energy-------
1941 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1942 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1943 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1944 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1945 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1946 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1947 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1948 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1952 !--------------------------------
1955 read (ifourier,*,end=115,err=115)
1956 read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13)
1958 write (iout,*) 'Type',i
1959 write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13)
1967 B1tilde(1,i) = buse(3)
1968 B1tilde(2,i) =-buse(5)
1969 B1tilde(1,-i) =-buse(3)
1970 B1tilde(2,-i) =buse(5)
1971 ! buse1tilde(1,i)=0.0d0
1972 ! buse1tilde(2,i)=0.0d0
1992 Ctilde(1,1,i)=buse(7)
1993 Ctilde(1,2,i)=buse(9)
1994 Ctilde(2,1,i)=-buse(9)
1995 Ctilde(2,2,i)=buse(7)
1996 Ctilde(1,1,-i)=buse(7)
1997 Ctilde(1,2,-i)=-buse(9)
1998 Ctilde(2,1,-i)=buse(9)
1999 Ctilde(2,2,-i)=buse(7)
2001 ! Ctilde(1,1,i)=0.0d0
2002 ! Ctilde(1,2,i)=0.0d0
2003 ! Ctilde(2,1,i)=0.0d0
2004 ! Ctilde(2,2,i)=0.0d0
2017 Dtilde(1,1,i)=buse(6)
2018 Dtilde(1,2,i)=buse(8)
2019 Dtilde(2,1,i)=-buse(8)
2020 Dtilde(2,2,i)=buse(6)
2021 Dtilde(1,1,-i)=buse(6)
2022 Dtilde(1,2,-i)=-buse(8)
2023 Dtilde(2,1,-i)=buse(8)
2024 Dtilde(2,2,-i)=buse(6)
2026 ! Dtilde(1,1,i)=0.0d0
2027 ! Dtilde(1,2,i)=0.0d0
2028 ! Dtilde(2,1,i)=0.0d0
2029 ! Dtilde(2,2,i)=0.0d0
2030 EE(1,1,i)= buse(10)+buse(11)
2031 EE(2,2,i)=-buse(10)+buse(11)
2032 EE(2,1,i)= buse(12)-buse(13)
2033 EE(1,2,i)= buse(12)+buse(13)
2034 EE(1,1,-i)= buse(10)+buse(11)
2035 EE(2,2,-i)=-buse(10)+buse(11)
2036 EE(2,1,-i)=-buse(12)+buse(13)
2037 EE(1,2,-i)=-buse(12)-buse(13)
2043 ! ee(2,1,i)=ee(1,2,i)
2047 write (iout,*) 'Type',i
2049 write(iout,*) B1(1,i),B1(2,i)
2051 write(iout,*) B2(1,i),B2(2,i)
2054 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2058 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2062 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2067 ! Read electrostatic-interaction parameters
2072 write (iout,'(/a)') 'Electrostatic interaction constants:'
2073 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2074 'IT','JT','APP','BPP','AEL6','AEL3'
2076 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2077 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2078 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2079 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2084 app (i,j)=epp(i,j)*rri*rri
2085 bpp (i,j)=-2.0D0*epp(i,j)*rri
2086 ael6(i,j)=elpp6(i,j)*4.2D0**6
2087 ael3(i,j)=elpp3(i,j)*4.2D0**3
2089 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2095 ! Read side-chain interaction parameters.
2097 !el from module energy - COMMON.INTERACT-------
2098 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2099 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2100 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2102 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2103 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2104 allocate(epslip(ntyp,ntyp))
2112 !--------------------------------
2114 read (isidep,*,end=117,err=117) ipot,expon
2115 if (ipot.lt.1 .or. ipot.gt.5) then
2116 write (iout,'(2a)') 'Error while reading SC interaction',&
2117 'potential file - unknown potential type.'
2119 call MPI_Finalize(Ierror)
2125 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2126 ', exponents are ',expon,2*expon
2127 ! goto (10,20,30,30,40) ipot
2129 !----------------------- LJ potential ---------------------------------
2131 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2132 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2133 (sigma0(i),i=1,ntyp)
2135 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2136 write (iout,'(a/)') 'The epsilon array:'
2137 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2138 write (iout,'(/a)') 'One-body parameters:'
2139 write (iout,'(a,4x,a)') 'residue','sigma'
2140 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2143 !----------------------- LJK potential --------------------------------
2145 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2146 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2147 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2149 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2150 write (iout,'(a/)') 'The epsilon array:'
2151 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2152 write (iout,'(/a)') 'One-body parameters:'
2153 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2154 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2158 !---------------------- GB or BP potential -----------------------------
2161 ! print *,"I AM in SCELE",scelemode
2162 if (scelemode.eq.0) then
2164 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2166 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2167 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2168 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2169 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2171 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2174 ! For the GB potential convert sigma'**2 into chi'
2177 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2181 write (iout,'(/a/)') 'Parameters of the BP potential:'
2182 write (iout,'(a/)') 'The epsilon array:'
2183 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2184 write (iout,'(/a)') 'One-body parameters:'
2185 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2187 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2188 chip(i),alp(i),i=1,ntyp)
2191 ! print *,ntyp,"NTYP"
2192 allocate(icharge(ntyp1))
2193 ! print *,ntyp,icharge(i)
2195 read (isidep,*) (icharge(i),i=1,ntyp)
2196 print *,ntyp,icharge(i)
2197 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2198 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2199 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2200 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2201 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2202 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2203 allocate(epsintab(ntyp,ntyp))
2204 allocate(dtail(2,ntyp,ntyp))
2205 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2206 allocate(wqdip(2,ntyp,ntyp))
2207 allocate(wstate(4,ntyp,ntyp))
2208 allocate(dhead(2,2,ntyp,ntyp))
2209 allocate(nstate(ntyp,ntyp))
2210 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2211 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2214 ! write (*,*) "Im in ALAB", i, " ", j
2216 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2217 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2218 chis(i,j),chis(j,i), &
2219 nstate(i,j),(wstate(k,i,j),k=1,4), &
2220 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2221 dtail(1,i,j),dtail(2,i,j), &
2222 epshead(i,j),sig0head(i,j), &
2223 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2224 alphapol(i,j),alphapol(j,i), &
2225 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
2226 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2232 sigma(i,j) = sigma(j,i)
2233 sigmap1(i,j)=sigmap1(j,i)
2234 sigmap2(i,j)=sigmap2(j,i)
2235 sigiso1(i,j)=sigiso1(j,i)
2236 sigiso2(i,j)=sigiso2(j,i)
2237 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2238 nstate(i,j) = nstate(j,i)
2239 dtail(1,i,j) = dtail(1,j,i)
2240 dtail(2,i,j) = dtail(2,j,i)
2242 alphasur(k,i,j) = alphasur(k,j,i)
2243 wstate(k,i,j) = wstate(k,j,i)
2244 alphiso(k,i,j) = alphiso(k,j,i)
2247 dhead(2,1,i,j) = dhead(1,1,j,i)
2248 dhead(2,2,i,j) = dhead(1,2,j,i)
2249 dhead(1,1,i,j) = dhead(2,1,j,i)
2250 dhead(1,2,i,j) = dhead(2,2,j,i)
2252 epshead(i,j) = epshead(j,i)
2253 sig0head(i,j) = sig0head(j,i)
2256 wqdip(k,i,j) = wqdip(k,j,i)
2259 wquad(i,j) = wquad(j,i)
2260 epsintab(i,j) = epsintab(j,i)
2261 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2266 !--------------------- GBV potential -----------------------------------
2268 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2269 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2270 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2271 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2273 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2274 write (iout,'(a/)') 'The epsilon array:'
2275 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2276 write (iout,'(/a)') 'One-body parameters:'
2277 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2278 's||/s_|_^2',' chip ',' alph '
2279 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2280 sigii(i),chip(i),alp(i),i=1,ntyp)
2283 write(iout,*)"Wrong ipot"
2289 !-----------------------------------------------------------------------
2290 ! Calculate the "working" parameters of SC interactions.
2292 !el from module energy - COMMON.INTERACT-------
2293 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2294 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2295 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2296 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2297 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2298 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2300 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2306 if (scelemode.eq.0) then
2315 sc_aa_tube_par(:)=0.0d0
2316 sc_bb_tube_par(:)=0.0d0
2318 !--------------------------------
2323 epslip(i,j)=epslip(j,i)
2326 if (scelemode.eq.0) then
2329 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2330 sigma(j,i)=sigma(i,j)
2331 rs0(i,j)=dwa16*sigma(i,j)
2336 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2337 'Working parameters of the SC interactions:',&
2338 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2343 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2345 ! print *,"SIGMA ZLA?",sigma(i,j)
2353 sigeps=dsign(1.0D0,epsij)
2355 aa_aq(i,j)=epsij*rrij*rrij
2356 print *,"ADASKO",epsij,rrij,expon
2357 bb_aq(i,j)=-sigeps*epsij*rrij
2358 aa_aq(j,i)=aa_aq(i,j)
2359 bb_aq(j,i)=bb_aq(i,j)
2360 epsijlip=epslip(i,j)
2361 sigeps=dsign(1.0D0,epsijlip)
2362 epsijlip=dabs(epsijlip)
2363 aa_lip(i,j)=epsijlip*rrij*rrij
2364 bb_lip(i,j)=-sigeps*epsijlip*rrij
2365 aa_lip(j,i)=aa_lip(i,j)
2366 bb_lip(j,i)=bb_lip(i,j)
2367 !C write(iout,*) aa_lip
2368 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2369 sigt1sq=sigma0(i)**2
2370 sigt2sq=sigma0(j)**2
2373 ratsig1=sigt2sq/sigt1sq
2374 ratsig2=1.0D0/ratsig1
2375 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2376 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2377 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2381 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2382 sigmaii(i,j)=rsum_max
2383 sigmaii(j,i)=rsum_max
2385 ! sigmaii(i,j)=r0(i,j)
2386 ! sigmaii(j,i)=r0(i,j)
2388 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2389 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2390 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2391 augm(i,j)=epsij*r_augm**(2*expon)
2392 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2399 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2400 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2401 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2406 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2407 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2408 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2409 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2410 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2411 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2412 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2413 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2414 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2415 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2416 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2417 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2418 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2419 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2420 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2421 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2430 read (isidep_nucl,*) ipot_nucl
2431 ! print *,"TU?!",ipot_nucl
2432 if (ipot_nucl.eq.1) then
2433 do i=1,ntyp_molec(2)
2434 do j=i,ntyp_molec(2)
2435 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2436 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2440 do i=1,ntyp_molec(2)
2441 do j=i,ntyp_molec(2)
2442 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2443 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2444 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2448 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2449 do i=1,ntyp_molec(2)
2450 do j=i,ntyp_molec(2)
2451 rrij=sigma_nucl(i,j)
2455 epsij=4*eps_nucl(i,j)
2456 sigeps=dsign(1.0D0,epsij)
2458 aa_nucl(i,j)=epsij*rrij*rrij
2459 bb_nucl(i,j)=-sigeps*epsij*rrij
2460 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2461 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2462 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2463 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2464 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2465 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2468 aa_nucl(i,j)=aa_nucl(j,i)
2469 bb_nucl(i,j)=bb_nucl(j,i)
2470 ael3_nucl(i,j)=ael3_nucl(j,i)
2471 ael6_nucl(i,j)=ael6_nucl(j,i)
2472 ael63_nucl(i,j)=ael63_nucl(j,i)
2473 ael32_nucl(i,j)=ael32_nucl(j,i)
2474 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2475 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2476 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2477 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2478 eps_nucl(i,j)=eps_nucl(j,i)
2479 sigma_nucl(i,j)=sigma_nucl(j,i)
2480 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2484 write(iout,*) "tube param"
2485 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2486 ccavtubpep,dcavtubpep,tubetranenepep
2487 sigmapeptube=sigmapeptube**6
2488 sigeps=dsign(1.0D0,epspeptube)
2489 epspeptube=dabs(epspeptube)
2490 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2491 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2492 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2494 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2495 ccavtub(i),dcavtub(i),tubetranene(i)
2496 sigmasctube=sigmasctube**6
2497 sigeps=dsign(1.0D0,epssctube)
2498 epssctube=dabs(epssctube)
2499 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2500 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2501 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2503 !-----------------READING SC BASE POTENTIALS-----------------------------
2504 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2505 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2506 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2507 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2508 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2509 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2510 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2511 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2512 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2513 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2514 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2515 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2516 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2517 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2518 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2519 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2522 do i=1,ntyp_molec(1)
2523 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2524 write (*,*) "Im in ", i, " ", j
2525 read(isidep_scbase,*) &
2526 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2527 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2528 write(*,*) "eps",eps_scbase(i,j)
2529 read(isidep_scbase,*) &
2530 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2531 chis_scbase(i,j,1),chis_scbase(i,j,2)
2532 read(isidep_scbase,*) &
2533 dhead_scbasei(i,j), &
2534 dhead_scbasej(i,j), &
2535 rborn_scbasei(i,j),rborn_scbasej(i,j)
2536 read(isidep_scbase,*) &
2537 (wdipdip_scbase(k,i,j),k=1,3), &
2538 (wqdip_scbase(k,i,j),k=1,2)
2539 read(isidep_scbase,*) &
2540 alphapol_scbase(i,j), &
2541 epsintab_scbase(i,j)
2544 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2545 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2547 do i=1,ntyp_molec(1)
2548 do j=1,ntyp_molec(2)-1
2549 epsij=eps_scbase(i,j)
2550 rrij=sigma_scbase(i,j)
2555 sigeps=dsign(1.0D0,epsij)
2557 aa_scbase(i,j)=epsij*rrij*rrij
2558 bb_scbase(i,j)=-sigeps*epsij*rrij
2561 !-----------------READING PEP BASE POTENTIALS-------------------
2562 allocate(eps_pepbase(ntyp_molec(2)))
2563 allocate(sigma_pepbase(ntyp_molec(2)))
2564 allocate(chi_pepbase(ntyp_molec(2),2))
2565 allocate(chipp_pepbase(ntyp_molec(2),2))
2566 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2567 allocate(sigmap1_pepbase(ntyp_molec(2)))
2568 allocate(sigmap2_pepbase(ntyp_molec(2)))
2569 allocate(chis_pepbase(ntyp_molec(2),2))
2570 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2573 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2574 write (*,*) "Im in ", i, " ", j
2575 read(isidep_pepbase,*) &
2576 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2577 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2578 write(*,*) "eps",eps_pepbase(j)
2579 read(isidep_pepbase,*) &
2580 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2581 chis_pepbase(j,1),chis_pepbase(j,2)
2582 read(isidep_pepbase,*) &
2583 (wdipdip_pepbase(k,j),k=1,3)
2585 allocate(aa_pepbase(ntyp_molec(2)))
2586 allocate(bb_pepbase(ntyp_molec(2)))
2588 do j=1,ntyp_molec(2)-1
2589 epsij=eps_pepbase(j)
2590 rrij=sigma_pepbase(j)
2595 sigeps=dsign(1.0D0,epsij)
2597 aa_pepbase(j)=epsij*rrij*rrij
2598 bb_pepbase(j)=-sigeps*epsij*rrij
2600 !--------------READING SC PHOSPHATE-------------------------------------
2601 allocate(eps_scpho(ntyp_molec(1)))
2602 allocate(sigma_scpho(ntyp_molec(1)))
2603 allocate(chi_scpho(ntyp_molec(1),2))
2604 allocate(chipp_scpho(ntyp_molec(1),2))
2605 allocate(alphasur_scpho(4,ntyp_molec(1)))
2606 allocate(sigmap1_scpho(ntyp_molec(1)))
2607 allocate(sigmap2_scpho(ntyp_molec(1)))
2608 allocate(chis_scpho(ntyp_molec(1),2))
2609 allocate(wqq_scpho(ntyp_molec(1)))
2610 allocate(wqdip_scpho(2,ntyp_molec(1)))
2611 allocate(alphapol_scpho(ntyp_molec(1)))
2612 allocate(epsintab_scpho(ntyp_molec(1)))
2613 allocate(dhead_scphoi(ntyp_molec(1)))
2614 allocate(rborn_scphoi(ntyp_molec(1)))
2615 allocate(rborn_scphoj(ntyp_molec(1)))
2616 allocate(alphi_scpho(ntyp_molec(1)))
2620 do j=1,ntyp_molec(1) ! without U then we will take T for U
2621 write (*,*) "Im in scpho ", i, " ", j
2622 read(isidep_scpho,*) &
2623 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2624 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2625 write(*,*) "eps",eps_scpho(j)
2626 read(isidep_scpho,*) &
2627 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2628 chis_scpho(j,1),chis_scpho(j,2)
2629 read(isidep_scpho,*) &
2630 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2631 read(isidep_scpho,*) &
2632 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2636 allocate(aa_scpho(ntyp_molec(1)))
2637 allocate(bb_scpho(ntyp_molec(1)))
2639 do j=1,ntyp_molec(1)
2646 sigeps=dsign(1.0D0,epsij)
2648 aa_scpho(j)=epsij*rrij*rrij
2649 bb_scpho(j)=-sigeps*epsij*rrij
2653 read(isidep_peppho,*) &
2654 eps_peppho,sigma_peppho
2655 read(isidep_peppho,*) &
2656 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2657 read(isidep_peppho,*) &
2658 (wqdip_peppho(k),k=1,2)
2666 sigeps=dsign(1.0D0,epsij)
2668 aa_peppho=epsij*rrij*rrij
2669 bb_peppho=-sigeps*epsij*rrij
2672 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2677 ! Define the SC-p interaction constants (hard-coded; old style)
2680 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2682 ! aad(i,1)=0.3D0*4.0D0**12
2683 ! Following line for constants currently implemented
2684 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2685 aad(i,1)=1.5D0*4.0D0**12
2686 ! aad(i,1)=0.17D0*5.6D0**12
2688 ! "Soft" SC-p repulsion
2690 ! Following line for constants currently implemented
2691 ! aad(i,1)=0.3D0*4.0D0**6
2692 ! "Hard" SC-p repulsion
2693 bad(i,1)=3.0D0*4.0D0**6
2694 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2703 ! 8/9/01 Read the SC-p interaction constants from file
2706 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2709 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2710 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2711 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2712 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2716 write (iout,*) "Parameters of SC-p interactions:"
2718 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2719 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2724 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2726 do i=1,ntyp_molec(2)
2727 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2729 do i=1,ntyp_molec(2)
2730 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2731 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2733 r0pp=1.12246204830937298142*5.16158
2739 ! Define the constants of the disulfide bridge
2743 ! Old arbitrary potential - commented out.
2748 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2749 ! energy surface of diethyl disulfide.
2750 ! A. Liwo and U. Kozlowska, 11/24/03
2767 write (iout,'(/a)') "Disulfide bridge parameters:"
2768 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2769 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2770 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2771 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2775 111 write (iout,*) "Error reading bending energy parameters."
2777 112 write (iout,*) "Error reading rotamer energy parameters."
2779 113 write (iout,*) "Error reading torsional energy parameters."
2781 114 write (iout,*) "Error reading double torsional energy parameters."
2783 115 write (iout,*) &
2784 "Error reading cumulant (multibody energy) parameters."
2786 116 write (iout,*) "Error reading electrostatic energy parameters."
2788 117 write (iout,*) "Error reading side chain interaction parameters."
2790 118 write (iout,*) "Error reading SCp interaction parameters."
2792 119 write (iout,*) "Error reading SCCOR parameters"
2795 call MPI_Finalize(Ierror)
2799 end subroutine parmread
2801 !-----------------------------------------------------------------------------
2803 !-----------------------------------------------------------------------------
2804 subroutine printmat(ldim,m,n,iout,key,a)
2807 character(len=3),dimension(n) :: key
2808 real(kind=8),dimension(ldim,n) :: a
2810 integer :: i,j,k,m,iout,nlim
2814 write (iout,1000) (key(k),k=i,nlim)
2816 1000 format (/5x,8(6x,a3))
2817 1020 format (/80(1h-)/)
2819 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2822 1010 format (a3,2x,8(f9.4))
2824 end subroutine printmat
2825 !-----------------------------------------------------------------------------
2827 !-----------------------------------------------------------------------------
2829 ! Read the PDB file and convert the peptide geometry into virtual-chain
2832 use energy_data, only: itype,istype
2836 ! use control, only: rescode,sugarcode
2837 ! implicit real*8 (a-h,o-z)
2838 ! include 'DIMENSIONS'
2839 ! include 'COMMON.LOCAL'
2840 ! include 'COMMON.VAR'
2841 ! include 'COMMON.CHAIN'
2842 ! include 'COMMON.INTERACT'
2843 ! include 'COMMON.IOUNITS'
2844 ! include 'COMMON.GEO'
2845 ! include 'COMMON.NAMES'
2846 ! include 'COMMON.CONTROL'
2847 ! include 'COMMON.DISTFIT'
2848 ! include 'COMMON.SETUP'
2849 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2851 logical :: lprn=.true.,fail
2852 real(kind=8),dimension(3) :: e1,e2,e3
2853 real(kind=8) :: dcj,efree_temp
2854 character(len=3) :: seq,res
2855 character(len=5) :: atom
2856 character(len=80) :: card
2857 real(kind=8),dimension(3,20) :: sccor
2858 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
2859 integer :: isugar,molecprev,firstion
2860 character*1 :: sugar
2862 real(kind=8),dimension(3) :: ccc
2864 integer,dimension(2,maxres/3) :: hfrag_alloc
2865 integer,dimension(4,maxres/3) :: bfrag_alloc
2866 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2867 real(kind=8),dimension(:,:), allocatable :: c_temporary
2868 integer,dimension(:,:) , allocatable :: itype_temporary
2869 integer,dimension(:),allocatable :: istype_temp
2876 ! write (2,*) "UNRES_PDB",unres_pdb
2896 !-----------------------------
2897 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2898 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2899 if(.not. allocated(istype)) allocate(istype(maxres))
2901 read (ipdbin,'(a80)',end=10) card
2902 write (iout,'(a)') card
2903 if (card(:5).eq.'HELIX') then
2906 read(card(22:25),*) hfrag(1,nhfrag)
2907 read(card(34:37),*) hfrag(2,nhfrag)
2909 if (card(:5).eq.'SHEET') then
2912 read(card(24:26),*) bfrag(1,nbfrag)
2913 read(card(35:37),*) bfrag(2,nbfrag)
2914 !rc----------------------------------------
2915 !rc to be corrected !!!
2916 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2917 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2918 !rc----------------------------------------
2920 if (card(:3).eq.'END') then
2922 else if (card(:3).eq.'TER') then
2927 itype(ires_old,molecule)=ntyp1_molec(molecule)
2928 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2929 nres_molec(molecule)=nres_molec(molecule)+2
2931 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2934 dc(j,ires)=sccor(j,iii)
2937 call sccenter(ires,iii,sccor)
2943 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2944 ! Fish out the ATOM cards.
2945 ! write(iout,*) 'card',card(1:20)
2946 if (index(card(1:4),'ATOM').gt.0) then
2947 read (card(12:16),*) atom
2948 ! write (iout,*) "! ",atom," !",ires
2949 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2950 read (card(23:26),*) ires
2951 read (card(18:20),'(a3)') res
2952 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2953 ! & " ires_old",ires_old
2954 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2955 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2956 if (ires-ishift+ishift1.ne.ires_old) then
2957 ! Calculate the CM of the preceding residue.
2958 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2960 ! write (iout,*) "Calculating sidechain center iii",iii
2963 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2966 call sccenter(ires_old,iii,sccor)
2970 ! Start new residue.
2971 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2974 else if (ibeg.eq.1) then
2975 write (iout,*) "BEG ires",ires
2977 if (res.ne.'GLY' .and. res.ne. 'ACE') then
2980 nres_molec(molecule)=nres_molec(molecule)+1
2982 ires=ires-ishift+ishift1
2984 ! write (iout,*) "ishift",ishift," ires",ires,&
2985 ! " ires_old",ires_old
2987 else if (ibeg.eq.2) then
2989 ishift=-ires_old+ires-1 !!!!!
2990 ishift1=ishift1-1 !!!!!
2991 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
2992 ires=ires-ishift+ishift1
2993 ! print *,ires,ishift,ishift1
2997 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
2998 ires=ires-ishift+ishift1
3001 ! print *,'atom',ires,atom
3002 if (res.eq.'ACE' .or. res.eq.'NHE') then
3005 if (atom.eq.'CA '.or.atom.eq.'N ') then
3007 itype(ires,molecule)=rescode(ires,res,0,molecule)
3009 ! nres_molec(molecule)=nres_molec(molecule)+1
3012 itype(ires,molecule)=rescode(ires,res(2:3),0,molecule)
3013 ! nres_molec(molecule)=nres_molec(molecule)+1
3014 read (card(19:19),'(a1)') sugar
3015 isugar=sugarcode(sugar,ires)
3016 ! if (ibeg.eq.1) then
3020 ! print *,"ires=",ires,istype(ires)
3026 ires=ires-ishift+ishift1
3028 ! write (iout,*) "ires_old",ires_old," ires",ires
3029 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3032 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3033 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3034 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3035 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3036 ! print *,ires,ishift,ishift1
3037 ! write (iout,*) "backbone ",atom
3039 write (iout,'(2i3,2x,a,3f8.3)') &
3040 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3043 nres_molec(molecule)=nres_molec(molecule)+1
3045 sccor(j,iii)=c(j,ires)
3047 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3048 atom.eq."C2'" .or. atom.eq."C3'" &
3049 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3050 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3051 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3052 ! print *,ires,ishift,ishift1
3056 ! sccor(j,iii)=c(j,ires)
3059 c(j,ires)=c(j,ires)+ccc(j)/5.0
3061 print *,counter,molecule
3062 if (counter.eq.5) then
3064 nres_molec(molecule)=nres_molec(molecule)+1
3067 ! sccor(j,iii)=c(j,ires)
3071 ! print *, "ATOM",atom(1:3)
3072 else if (atom.eq."C5'") then
3073 read (card(19:19),'(a1)') sugar
3074 isugar=sugarcode(sugar,ires)
3079 ! print *,ires,istype(ires)
3083 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3084 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3085 nres_molec(molecule)=nres_molec(molecule)+1
3086 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3090 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3092 else if ((atom.eq."C1'").and.unres_pdb) then
3094 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3095 ! write (*,*) card(23:27),ires,itype(ires,1)
3096 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3097 atom.ne.'N' .and. atom.ne.'C' .and. &
3098 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3099 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3100 .and. atom.ne.'P '.and. &
3101 atom(1:1).ne.'H' .and. &
3102 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3104 ! write (iout,*) "sidechain ",atom
3105 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3106 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3107 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3109 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3112 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3113 if (firstion.eq.0) then
3117 dc(j,ires)=sccor(j,iii)
3120 call sccenter(ires,iii,sccor)
3123 read (card(12:16),*) atom
3124 print *,"HETATOM", atom
3125 read (card(18:20),'(a3)') res
3126 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3127 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3128 .or.(atom(1:2).eq.'K ')) &
3131 if (molecule.ne.5) molecprev=molecule
3133 nres_molec(molecule)=nres_molec(molecule)+1
3134 print *,"HERE",nres_molec(molecule)
3136 itype(ires,molecule)=rescode(ires,res,0,molecule)
3137 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3141 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3142 if (ires.eq.0) return
3143 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3146 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3148 nres_molec(molecule)=nres_molec(molecule)-2
3149 print *,'I have',nres, nres_molec(:)
3151 do k=1,4 ! ions are without dummy
3152 if (nres_molec(k).eq.0) cycle
3154 ! write (iout,*) i,itype(i,1)
3155 ! if (itype(i,1).eq.ntyp1) then
3156 ! write (iout,*) "dummy",i,itype(i,1)
3158 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3159 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3163 if (itype(i,k).eq.ntyp1_molec(k)) then
3164 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3165 if (itype(i+2,k).eq.0) then
3166 ! print *,"masz sieczke"
3168 if (itype(i+2,j).ne.0) then
3170 itype(i+1,j)=ntyp1_molec(j)
3171 nres_molec(k)=nres_molec(k)-1
3172 nres_molec(j)=nres_molec(j)+1
3178 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3179 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3180 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3181 ! if (unres_pdb) then
3182 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3183 ! print *,i,'tu dochodze'
3184 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3192 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3196 dcj=(c(j,i-2)-c(j,i-3))/2.0
3197 if (dcj.eq.0) dcj=1.23591524223
3202 else !itype(i+1,1).eq.ntyp1
3203 ! if (unres_pdb) then
3204 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3205 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3212 c(j,i)=c(j,i+1)-1.9d0*e2(j)
3216 dcj=(c(j,i+3)-c(j,i+2))/2.0
3217 if (dcj.eq.0) dcj=1.23591524223
3222 endif !itype(i+1,1).eq.ntyp1
3223 endif !itype.eq.ntyp1
3227 ! Calculate the CM of the last side chain.
3231 dc(j,ires)=sccor(j,iii)
3234 call sccenter(ires,iii,sccor)
3240 ! print *,"molecule",molecule
3241 if ((itype(nres,1).ne.10)) then
3243 if (molecule.eq.5) molecule=molecprev
3244 itype(nres,molecule)=ntyp1_molec(molecule)
3245 nres_molec(molecule)=nres_molec(molecule)+1
3246 ! if (unres_pdb) then
3247 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3248 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3255 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3259 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3260 c(j,nres)=c(j,nres-1)+dcj
3261 c(j,2*nres)=c(j,nres)
3265 ! print *,'I have',nres, nres_molec(:)
3267 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3269 if (nres.ne.nres0) then
3270 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3272 stop "Error nres value in WHAM input"
3275 !---------------------------------
3276 !el reallocate tables
3279 ! hfrag_alloc(j,i)=hfrag(j,i)
3282 ! bfrag_alloc(j,i)=bfrag(j,i)
3288 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3289 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3290 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3291 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3295 ! hfrag(j,i)=hfrag_alloc(j,i)
3300 ! bfrag(j,i)=bfrag_alloc(j,i)
3303 !el end reallocate tables
3304 !---------------------------------
3312 c(j,2*nres)=c(j,nres)
3315 if (itype(1,1).eq.ntyp1) then
3319 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3320 call refsys(2,3,4,e1,e2,e3,fail)
3327 c(j,1)=c(j,2)-1.9d0*e2(j)
3331 dcj=(c(j,4)-c(j,3))/2.0
3337 ! First lets assign correct dummy to correct type of chain
3339 if (itype(1,1).eq.ntyp1) then
3340 if (itype(2,1).eq.0) then
3342 if (itype(2,j).ne.0) then
3344 itype(1,j)=ntyp1_molec(j)
3345 nres_molec(1)=nres_molec(1)-1
3346 nres_molec(j)=nres_molec(j)+1
3353 print *,'I have',nres, nres_molec(:)
3355 ! Copy the coordinates to reference coordinates
3361 ! Calculate internal coordinates.
3363 write (iout,'(/a)') &
3364 "Cartesian coordinates of the reference structure"
3365 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3366 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3368 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3369 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3370 (c(j,ires+nres),j=1,3)
3373 ! znamy już nres wiec mozna alokowac tablice
3374 ! Calculate internal coordinates.
3375 if(me.eq.king.or..not.out1file)then
3376 write (iout,'(a)') &
3377 "Backbone and SC coordinates as read from the PDB"
3379 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
3380 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3381 (c(j,nres+ires),j=1,3)
3384 ! NOW LETS ROCK! SORTING
3385 allocate(c_temporary(3,2*nres))
3386 allocate(itype_temporary(nres,5))
3387 if (.not.allocated(molnum)) allocate(molnum(nres+1))
3388 if (.not.allocated(istype)) write(iout,*) &
3389 "SOMETHING WRONG WITH ISTYTPE"
3390 allocate(istype_temp(nres))
3391 itype_temporary(:,:)=0
3395 if (itype(i,k).ne.0) then
3397 c_temporary(j,seqalingbegin)=c(j,i)
3398 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3401 itype_temporary(seqalingbegin,k)=itype(i,k)
3402 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3403 istype_temp(seqalingbegin)=istype(i)
3404 molnum(seqalingbegin)=k
3405 seqalingbegin=seqalingbegin+1
3411 c(j,i)=c_temporary(j,i)
3416 itype(i,k)=itype_temporary(i,k)
3417 istype(i)=istype_temp(i)
3420 ! if (itype(1,1).eq.ntyp1) then
3423 ! if (unres_pdb) then
3424 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3425 ! call refsys(2,3,4,e1,e2,e3,fail)
3432 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3436 ! dcj=(c(j,4)-c(j,3))/2.0
3438 ! c(j,nres+1)=c(j,1)
3444 write (iout,'(/a)') &
3445 "Cartesian coordinates of the reference structure after sorting"
3446 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3447 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3449 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3450 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3451 (c(j,ires+nres),j=1,3)
3455 ! print *,seqalingbegin,nres
3456 if(.not.allocated(vbld)) then
3457 allocate(vbld(2*nres))
3462 if(.not.allocated(vbld_inv)) then
3463 allocate(vbld_inv(2*nres))
3469 if(.not.allocated(theta)) then
3470 allocate(theta(nres+2))
3474 if(.not.allocated(phi)) allocate(phi(nres+2))
3475 if(.not.allocated(alph)) allocate(alph(nres+2))
3476 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3477 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3478 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3479 if(.not.allocated(costtab)) allocate(costtab(nres))
3480 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3481 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3482 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3483 if(.not.allocated(xxref)) allocate(xxref(nres))
3484 if(.not.allocated(yyref)) allocate(yyref(nres))
3485 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3486 if(.not.allocated(dc_norm)) then
3487 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3488 allocate(dc_norm(3,0:2*nres+2))
3492 call int_from_cart(.true.,.false.)
3493 call sc_loc_geom(.false.)
3495 thetaref(i)=theta(i)
3505 dc(j,i)=c(j,i+1)-c(j,i)
3506 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3511 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3512 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3514 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3518 ! Copy the coordinates to reference coordinates
3519 ! Splits to single chain if occurs
3523 ! cref(j,i,cou)=c(j,i)
3527 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3528 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3529 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3530 !-----------------------------
3534 write (iout,*) "symetr", symetr
3537 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3539 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3542 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3547 cref(j,i,cou)=c(j,i)
3548 cref(j,i+nres,cou)=c(j,i+nres)
3550 chain_rep(j,lll,kkk)=c(j,i)
3551 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3555 write (iout,*) chain_length
3556 if (chain_length.eq.0) chain_length=nres
3558 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3559 chain_rep(j,chain_length+nres,symetr) &
3560 =chain_rep(j,chain_length+nres,1)
3563 ! write (iout,*) "spraw lancuchy",chain_length,symetr
3565 ! do kkk=1,chain_length
3566 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
3570 ! makes copy of chains
3571 write (iout,*) "symetr", symetr
3576 if (symetr.gt.1) then
3583 write(iout,*) (tabperm(i,kkk),kkk=1,4)
3589 write (iout,*) i,icha
3590 do lll=1,chain_length
3592 if (cou.le.nres) then
3594 kupa=mod(lll,chain_length)
3595 iprzes=(kkk-1)*chain_length+lll
3596 if (kupa.eq.0) kupa=chain_length
3597 write (iout,*) "kupa", kupa
3598 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3599 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3606 !-koniec robienia kopii
3609 write (iout,*) "nowa struktura", nperm
3611 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3613 cref(3,i,kkk),cref(1,nres+i,kkk),&
3614 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3616 100 format (//' alpha-carbon coordinates ',&
3617 ' centroid coordinates'/ &
3618 ' ', 6X,'X',11X,'Y',11X,'Z', &
3619 10X,'X',11X,'Y',11X,'Z')
3620 110 format (a,'(',i5,')',6f12.5)
3626 bfrag(i,j)=bfrag(i,j)-ishift
3632 hfrag(i,j)=hfrag(i,j)-ishift
3638 end subroutine readpdb
3639 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3640 !-----------------------------------------------------------------------------
3642 !-----------------------------------------------------------------------------
3643 subroutine read_control
3657 use random, only: random_init
3658 ! implicit real*8 (a-h,o-z)
3659 ! include 'DIMENSIONS'
3661 use prng, only:prng_restart
3663 logical :: OKRandom!, prng_restart
3666 ! include 'COMMON.IOUNITS'
3667 ! include 'COMMON.TIME1'
3668 ! include 'COMMON.THREAD'
3669 ! include 'COMMON.SBRIDGE'
3670 ! include 'COMMON.CONTROL'
3671 ! include 'COMMON.MCM'
3672 ! include 'COMMON.MAP'
3673 ! include 'COMMON.HEADER'
3674 ! include 'COMMON.CSA'
3675 ! include 'COMMON.CHAIN'
3676 ! include 'COMMON.MUCA'
3677 ! include 'COMMON.MD'
3678 ! include 'COMMON.FFIELD'
3679 ! include 'COMMON.INTERACT'
3680 ! include 'COMMON.SETUP'
3681 !el integer :: KDIAG,ICORFL,IXDR
3682 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3683 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3684 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3685 ! character(len=80) :: ucase
3686 character(len=640) :: controlcard
3688 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3694 read (INP,'(a)') titel
3695 call card_concat(controlcard,.true.)
3696 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3697 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3698 call reada(controlcard,'SEED',seed,0.0D0)
3699 call random_init(seed)
3700 ! Set up the time limit (caution! The time must be input in minutes!)
3701 read_cart=index(controlcard,'READ_CART').gt.0
3702 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3703 call readi(controlcard,'SYM',symetr,1)
3704 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3705 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3706 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3707 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3708 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3709 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3710 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3711 call reada(controlcard,'DRMS',drms,0.1D0)
3712 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3713 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3714 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3715 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3716 write (iout,'(a,f10.1)')'DRMS = ',drms
3717 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3718 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3720 call readi(controlcard,'NZ_START',nz_start,0)
3721 call readi(controlcard,'NZ_END',nz_end,0)
3722 call readi(controlcard,'IZ_SC',iz_sc,0)
3723 timlim=60.0D0*timlim
3724 safety = 60.0d0*safety
3727 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3728 !C SHIELD keyword sets if the shielding effect of side-chains is used
3729 !C 0 denots no shielding is used all peptide are equally despite the
3730 !C solvent accesible area
3731 !C 1 the newly introduced function
3732 !C 2 reseved for further possible developement
3733 call readi(controlcard,'SHIELD',shield_mode,0)
3734 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3735 write(iout,*) "shield_mode",shield_mode
3736 !C Varibles set size of box
3737 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3738 protein=index(controlcard,"PROTEIN").gt.0
3739 ions=index(controlcard,"IONS").gt.0
3740 nucleic=index(controlcard,"NUCLEIC").gt.0
3741 write (iout,*) "with_theta_constr ",with_theta_constr
3742 AFMlog=(index(controlcard,'AFM'))
3743 selfguide=(index(controlcard,'SELFGUIDE'))
3744 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3745 call readi(controlcard,'GENCONSTR',genconstr,0)
3746 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3747 call reada(controlcard,'BOXY',boxysize,100.0d0)
3748 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3749 call readi(controlcard,'TUBEMOD',tubemode,0)
3750 print *,"SCELE",scelemode
3751 call readi(controlcard,"SCELEMODE",scelemode,0)
3752 print *,"SCELE",scelemode
3754 ! elemode = 0 is orignal UNRES electrostatics
3755 ! elemode = 1 is "Momo" potentials in progress
3756 ! elemode = 2 is in development EVALD
3757 write (iout,*) TUBEmode,"TUBEMODE"
3758 if (TUBEmode.gt.0) then
3759 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3760 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3761 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3762 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3763 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3764 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3765 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3766 buftubebot=bordtubebot+tubebufthick
3767 buftubetop=bordtubetop-tubebufthick
3770 ! CUTOFFF ON ELECTROSTATICS
3771 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3772 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3773 write(iout,*) "R_CUT_ELE=",r_cut_ele
3774 ! Lipidic parameters
3775 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3776 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3777 if (lipthick.gt.0.0d0) then
3778 bordliptop=(boxzsize+lipthick)/2.0
3779 bordlipbot=bordliptop-lipthick
3780 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3781 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3782 buflipbot=bordlipbot+lipbufthick
3783 bufliptop=bordliptop-lipbufthick
3784 if ((lipbufthick*2.0d0).gt.lipthick) &
3785 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3786 endif !lipthick.gt.0
3787 write(iout,*) "bordliptop=",bordliptop
3788 write(iout,*) "bordlipbot=",bordlipbot
3789 write(iout,*) "bufliptop=",bufliptop
3790 write(iout,*) "buflipbot=",buflipbot
3791 write (iout,*) "SHIELD MODE",shield_mode
3793 !C-------------------------
3794 minim=(index(controlcard,'MINIMIZE').gt.0)
3795 dccart=(index(controlcard,'CART').gt.0)
3796 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3797 overlapsc=.not.overlapsc
3798 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3799 searchsc=.not.searchsc
3800 sideadd=(index(controlcard,'SIDEADD').gt.0)
3801 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3802 outpdb=(index(controlcard,'PDBOUT').gt.0)
3803 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3804 pdbref=(index(controlcard,'PDBREF').gt.0)
3805 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3806 indpdb=index(controlcard,'PDBSTART')
3807 extconf=(index(controlcard,'EXTCONF').gt.0)
3808 call readi(controlcard,'IPRINT',iprint,0)
3809 call readi(controlcard,'MAXGEN',maxgen,10000)
3810 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3811 call readi(controlcard,"KDIAG",kdiag,0)
3812 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3813 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3814 write (iout,*) "RESCALE_MODE",rescale_mode
3815 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3816 if (index(controlcard,'REGULAR').gt.0.0D0) then
3817 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3821 if (index(controlcard,'CHECKGRAD').gt.0) then
3823 if (index(controlcard,'CART').gt.0) then
3825 elseif (index(controlcard,'CARINT').gt.0) then
3830 elseif (index(controlcard,'THREAD').gt.0) then
3832 call readi(controlcard,'THREAD',nthread,0)
3833 if (nthread.gt.0) then
3834 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3837 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3838 stop 'Error termination in Read_Control.'
3840 else if (index(controlcard,'MCMA').gt.0) then
3842 else if (index(controlcard,'MCEE').gt.0) then
3844 else if (index(controlcard,'MULTCONF').gt.0) then
3846 else if (index(controlcard,'MAP').gt.0) then
3848 call readi(controlcard,'MAP',nmap,0)
3849 else if (index(controlcard,'CSA').gt.0) then
3851 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3853 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3856 !fcm else if (index(controlcard,'MCMF').gt.0) then
3858 else if (index(controlcard,'SOFTREG').gt.0) then
3860 else if (index(controlcard,'CHECK_BOND').gt.0) then
3862 else if (index(controlcard,'TEST').gt.0) then
3864 else if (index(controlcard,'MD').gt.0) then
3866 else if (index(controlcard,'RE ').gt.0) then
3870 lmuca=index(controlcard,'MUCA').gt.0
3871 call readi(controlcard,'MUCADYN',mucadyn,0)
3872 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3873 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3875 write (iout,*) 'MUCADYN=',mucadyn
3876 write (iout,*) 'MUCASMOOTH=',muca_smooth
3879 iscode=index(controlcard,'ONE_LETTER')
3880 indphi=index(controlcard,'PHI')
3881 indback=index(controlcard,'BACK')
3882 iranconf=index(controlcard,'RAND_CONF')
3883 i2ndstr=index(controlcard,'USE_SEC_PRED')
3884 gradout=index(controlcard,'GRADOUT').gt.0
3885 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3886 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3887 if (me.eq.king .or. .not.out1file ) &
3888 write (iout,*) "DISTCHAINMAX",distchainmax
3890 if(me.eq.king.or..not.out1file) &
3891 write (iout,'(2a)') diagmeth(kdiag),&
3892 ' routine used to diagonalize matrices.'
3893 if (shield_mode.gt.0) then
3895 !C VSolvSphere the volume of solving sphere
3897 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3898 !C there will be no distinction between proline peptide group and normal peptide
3899 !C group in case of shielding parameters
3900 VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
3901 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
3902 write (iout,*) VSolvSphere,VSolvSphere_div
3903 !C long axis of side chain
3905 long_r_sidechain(i)=vbldsc0(1,i)
3906 short_r_sidechain(i)=sigma0(i)
3907 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3913 end subroutine read_control
3914 !-----------------------------------------------------------------------------
3915 subroutine read_REMDpar
3917 ! Read REMD settings
3924 use control_data, only:out1file
3925 ! implicit real*8 (a-h,o-z)
3926 ! include 'DIMENSIONS'
3927 ! include 'COMMON.IOUNITS'
3928 ! include 'COMMON.TIME1'
3929 ! include 'COMMON.MD'
3932 !el include 'COMMON.LANGEVIN'
3934 !el include 'COMMON.LANGEVIN.lang0'
3936 ! include 'COMMON.INTERACT'
3937 ! include 'COMMON.NAMES'
3938 ! include 'COMMON.GEO'
3939 ! include 'COMMON.REMD'
3940 ! include 'COMMON.CONTROL'
3941 ! include 'COMMON.SETUP'
3942 ! character(len=80) :: ucase
3943 character(len=320) :: controlcard
3944 character(len=3200) :: controlcard1
3945 integer :: iremd_m_total
3948 ! real(kind=8) :: var,ene
3950 if(me.eq.king.or..not.out1file) &
3951 write (iout,*) "REMD setup"
3953 call card_concat(controlcard,.true.)
3954 call readi(controlcard,"NREP",nrep,3)
3955 call readi(controlcard,"NSTEX",nstex,1000)
3956 call reada(controlcard,"RETMIN",retmin,10.0d0)
3957 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3958 mremdsync=(index(controlcard,'SYNC').gt.0)
3959 call readi(controlcard,"NSYN",i_sync_step,100)
3960 restart1file=(index(controlcard,'REST1FILE').gt.0)
3961 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3962 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3963 if(max_cache_traj_use.gt.max_cache_traj) &
3964 max_cache_traj_use=max_cache_traj
3965 if(me.eq.king.or..not.out1file) then
3966 !d if (traj1file) then
3967 !rc caching is in testing - NTWX is not ignored
3968 !d write (iout,*) "NTWX value is ignored"
3969 !d write (iout,*) " trajectory is stored to one file by master"
3970 !d write (iout,*) " before exchange at NSTEX intervals"
3972 write (iout,*) "NREP= ",nrep
3973 write (iout,*) "NSTEX= ",nstex
3974 write (iout,*) "SYNC= ",mremdsync
3975 write (iout,*) "NSYN= ",i_sync_step
3976 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
3979 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
3980 if (index(controlcard,'TLIST').gt.0) then
3982 call card_concat(controlcard1,.true.)
3983 read(controlcard1,*) (remd_t(i),i=1,nrep)
3984 if(me.eq.king.or..not.out1file) &
3985 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
3988 if (index(controlcard,'MLIST').gt.0) then
3990 call card_concat(controlcard1,.true.)
3991 read(controlcard1,*) (remd_m(i),i=1,nrep)
3992 if(me.eq.king.or..not.out1file) then
3993 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
3996 iremd_m_total=iremd_m_total+remd_m(i)
3998 write (iout,*) 'Total number of replicas ',iremd_m_total
4001 if(me.eq.king.or..not.out1file) &
4002 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4004 end subroutine read_REMDpar
4005 !-----------------------------------------------------------------------------
4006 subroutine read_MDpar
4010 use control_data, only: r_cut,rlamb,out1file
4012 use geometry_data, only: pi
4014 ! implicit real*8 (a-h,o-z)
4015 ! include 'DIMENSIONS'
4016 ! include 'COMMON.IOUNITS'
4017 ! include 'COMMON.TIME1'
4018 ! include 'COMMON.MD'
4021 !el include 'COMMON.LANGEVIN'
4023 !el include 'COMMON.LANGEVIN.lang0'
4025 ! include 'COMMON.INTERACT'
4026 ! include 'COMMON.NAMES'
4027 ! include 'COMMON.GEO'
4028 ! include 'COMMON.SETUP'
4029 ! include 'COMMON.CONTROL'
4030 ! include 'COMMON.SPLITELE'
4031 ! character(len=80) :: ucase
4032 character(len=320) :: controlcard
4037 call card_concat(controlcard,.true.)
4038 call readi(controlcard,"NSTEP",n_timestep,1000000)
4039 call readi(controlcard,"NTWE",ntwe,100)
4040 call readi(controlcard,"NTWX",ntwx,1000)
4041 call reada(controlcard,"DT",d_time,1.0d-1)
4042 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4043 call reada(controlcard,"DAMAX",damax,1.0d1)
4044 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4045 call readi(controlcard,"LANG",lang,0)
4046 RESPA = index(controlcard,"RESPA") .gt. 0
4047 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4048 ntime_split0=ntime_split
4049 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4050 ntime_split0=ntime_split
4051 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4052 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4053 rest = index(controlcard,"REST").gt.0
4054 tbf = index(controlcard,"TBF").gt.0
4055 usampl = index(controlcard,"USAMPL").gt.0
4056 mdpdb = index(controlcard,"MDPDB").gt.0
4057 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4058 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4059 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4060 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4061 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4062 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4063 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4064 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4065 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4066 large = index(controlcard,"LARGE").gt.0
4067 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4068 rattle = index(controlcard,"RATTLE").gt.0
4069 preminim=(index(controlcard,'PREMINIM').gt.0)
4070 write (iout,*) "PREMINIM ",preminim
4071 dccart=(index(controlcard,'CART').gt.0)
4072 if (preminim) call read_minim
4073 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4079 if(me.eq.king.or..not.out1file) then
4081 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4083 write (iout,'(a)') "The units are:"
4084 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4085 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4086 " acceleration: angstrom/(48.9 fs)**2"
4087 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4089 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4090 write (iout,'(a60,f10.5,a)') &
4091 "Initial time step of numerical integration:",d_time,&
4093 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4095 write (iout,'(2a,i4,a)') &
4096 "A-MTS algorithm used; initial time step for fast-varying",&
4097 " short-range forces split into",ntime_split," steps."
4098 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4099 r_cut," lambda",rlamb
4101 write (iout,'(2a,f10.5)') &
4102 "Maximum acceleration threshold to reduce the time step",&
4103 "/increase split number:",damax
4104 write (iout,'(2a,f10.5)') &
4105 "Maximum predicted energy drift to reduce the timestep",&
4106 "/increase split number:",edriftmax
4107 write (iout,'(a60,f10.5)') &
4108 "Maximum velocity threshold to reduce velocities:",dvmax
4109 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4110 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4111 if (rattle) write (iout,'(a60)') &
4112 "Rattle algorithm used to constrain the virtual bonds"
4116 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4117 call reada(controlcard,"RWAT",rwat,1.4d0)
4118 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4119 surfarea=index(controlcard,"SURFAREA").gt.0
4120 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4121 if(me.eq.king.or..not.out1file)then
4122 write (iout,'(/a,$)') "Langevin dynamics calculation"
4124 write (iout,'(a/)') &
4125 " with direct integration of Langevin equations"
4126 else if (lang.eq.2) then
4127 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4128 else if (lang.eq.3) then
4129 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4130 else if (lang.eq.4) then
4131 write (iout,'(a/)') " in overdamped mode"
4133 write (iout,'(//a,i5)') &
4134 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4137 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4138 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4139 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4140 write (iout,'(a60,f10.5)') &
4141 "Scaling factor of the friction forces:",scal_fric
4142 if (surfarea) write (iout,'(2a,i10,a)') &
4143 "Friction coefficients will be scaled by solvent-accessible",&
4144 " surface area every",reset_fricmat," steps."
4146 ! Calculate friction coefficients and bounds of stochastic forces
4147 eta=6*pi*cPoise*etawat
4148 if(me.eq.king.or..not.out1file) &
4149 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4152 do j=1,5 !types of molecules
4153 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4154 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4156 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4157 do j=1,5 !types of molecules
4159 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4160 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4164 if(me.eq.king.or..not.out1file)then
4165 write (iout,'(/2a/)') &
4166 "Radii of site types and friction coefficients and std's of",&
4167 " stochastic forces of fully exposed sites"
4168 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4170 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4171 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4175 if(me.eq.king.or..not.out1file)then
4176 write (iout,'(a)') "Berendsen bath calculation"
4177 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4178 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4180 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4181 count_reset_moment," steps"
4183 write (iout,'(a,i10,a)') &
4184 "Velocities will be reset at random every",count_reset_vel,&
4188 if(me.eq.king.or..not.out1file) &
4189 write (iout,'(a31)') "Microcanonical mode calculation"
4191 if(me.eq.king.or..not.out1file)then
4192 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4194 write(iout,*) "MD running with constraints."
4195 write(iout,*) "Equilibration time ", eq_time, " mtus."
4196 write(iout,*) "Constraining ", nfrag," fragments."
4197 write(iout,*) "Length of each fragment, weight and q0:"
4199 write (iout,*) "Set of restraints #",iset
4201 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4202 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4204 write(iout,*) "constraints between ", npair, "fragments."
4205 write(iout,*) "constraint pairs, weights and q0:"
4207 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4208 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4210 write(iout,*) "angle constraints within ", nfrag_back,&
4211 "backbone fragments."
4212 write(iout,*) "fragment, weights:"
4214 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4215 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4216 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4219 iset=mod(kolor,nset)+1
4222 if(me.eq.king.or..not.out1file) &
4223 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4225 end subroutine read_MDpar
4226 !-----------------------------------------------------------------------------
4230 ! implicit real*8 (a-h,o-z)
4231 ! include 'DIMENSIONS'
4232 ! include 'COMMON.MAP'
4233 ! include 'COMMON.IOUNITS'
4234 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4235 character(len=80) :: mapcard !,ucase
4238 ! real(kind=8) :: var,ene
4241 read (inp,'(a)') mapcard
4242 mapcard=ucase(mapcard)
4243 if (index(mapcard,'PHI').gt.0) then
4245 else if (index(mapcard,'THE').gt.0) then
4247 else if (index(mapcard,'ALP').gt.0) then
4249 else if (index(mapcard,'OME').gt.0) then
4252 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4253 stop 'Error - illegal variable spec in MAP card.'
4255 call readi (mapcard,'RES1',res1(imap),0)
4256 call readi (mapcard,'RES2',res2(imap),0)
4257 if (res1(imap).eq.0) then
4258 res1(imap)=res2(imap)
4259 else if (res2(imap).eq.0) then
4260 res2(imap)=res1(imap)
4262 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4263 write (iout,'(a)') &
4264 'Error - illegal definition of variable group in MAP.'
4265 stop 'Error - illegal definition of variable group in MAP.'
4267 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4268 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4269 call readi(mapcard,'NSTEP',nstep(imap),0)
4270 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4271 write (iout,'(a)') &
4272 'Illegal boundary and/or step size specification in MAP.'
4273 stop 'Illegal boundary and/or step size specification in MAP.'
4277 end subroutine map_read
4278 !-----------------------------------------------------------------------------
4281 use control_data, only: vdisulf
4283 ! implicit real*8 (a-h,o-z)
4284 ! include 'DIMENSIONS'
4285 ! include 'COMMON.IOUNITS'
4286 ! include 'COMMON.GEO'
4287 ! include 'COMMON.CSA'
4288 ! include 'COMMON.BANK'
4289 ! include 'COMMON.CONTROL'
4290 ! character(len=80) :: ucase
4291 character(len=620) :: mcmcard
4293 ! integer :: ntf,ik,iw_pdb
4294 ! real(kind=8) :: var,ene
4296 call card_concat(mcmcard,.true.)
4298 call readi(mcmcard,'NCONF',nconf,50)
4299 call readi(mcmcard,'NADD',nadd,0)
4300 call readi(mcmcard,'JSTART',jstart,1)
4301 call readi(mcmcard,'JEND',jend,1)
4302 call readi(mcmcard,'NSTMAX',nstmax,500000)
4303 call readi(mcmcard,'N0',n0,1)
4304 call readi(mcmcard,'N1',n1,6)
4305 call readi(mcmcard,'N2',n2,4)
4306 call readi(mcmcard,'N3',n3,0)
4307 call readi(mcmcard,'N4',n4,0)
4308 call readi(mcmcard,'N5',n5,0)
4309 call readi(mcmcard,'N6',n6,10)
4310 call readi(mcmcard,'N7',n7,0)
4311 call readi(mcmcard,'N8',n8,0)
4312 call readi(mcmcard,'N9',n9,0)
4313 call readi(mcmcard,'N14',n14,0)
4314 call readi(mcmcard,'N15',n15,0)
4315 call readi(mcmcard,'N16',n16,0)
4316 call readi(mcmcard,'N17',n17,0)
4317 call readi(mcmcard,'N18',n18,0)
4319 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4321 call readi(mcmcard,'NDIFF',ndiff,2)
4322 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4323 call readi(mcmcard,'IS1',is1,1)
4324 call readi(mcmcard,'IS2',is2,8)
4325 call readi(mcmcard,'NRAN0',nran0,4)
4326 call readi(mcmcard,'NRAN1',nran1,2)
4327 call readi(mcmcard,'IRR',irr,1)
4328 call readi(mcmcard,'NSEED',nseed,20)
4329 call readi(mcmcard,'NTOTAL',ntotal,10000)
4330 call reada(mcmcard,'CUT1',cut1,2.0d0)
4331 call reada(mcmcard,'CUT2',cut2,5.0d0)
4332 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4333 call readi(mcmcard,'ICMAX',icmax,3)
4334 call readi(mcmcard,'IRESTART',irestart,0)
4335 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4338 call reada(mcmcard,'DELE',dele,20.0d0)
4339 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4340 call readi(mcmcard,'IREF',iref,0)
4341 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4342 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4343 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4344 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4345 write (iout,*) "NCONF_IN",nconf_in
4347 end subroutine csaread
4348 !-----------------------------------------------------------------------------
4352 use control_data, only: MaxMoveType
4355 ! implicit real*8 (a-h,o-z)
4356 ! include 'DIMENSIONS'
4357 ! include 'COMMON.MCM'
4358 ! include 'COMMON.MCE'
4359 ! include 'COMMON.IOUNITS'
4360 ! character(len=80) :: ucase
4361 character(len=320) :: mcmcard
4364 ! real(kind=8) :: var,ene
4366 call card_concat(mcmcard,.true.)
4367 call readi(mcmcard,'MAXACC',maxacc,100)
4368 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4369 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4370 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4371 call readi(mcmcard,'MAXREPM',maxrepm,200)
4372 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4373 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4374 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4375 call reada(mcmcard,'E_UP',e_up,5.0D0)
4376 call reada(mcmcard,'DELTE',delte,0.1D0)
4377 call readi(mcmcard,'NSWEEP',nsweep,5)
4378 call readi(mcmcard,'NSTEPH',nsteph,0)
4379 call readi(mcmcard,'NSTEPC',nstepc,0)
4380 call reada(mcmcard,'TMIN',tmin,298.0D0)
4381 call reada(mcmcard,'TMAX',tmax,298.0D0)
4382 call readi(mcmcard,'NWINDOW',nwindow,0)
4383 call readi(mcmcard,'PRINT_MC',print_mc,0)
4384 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4385 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4386 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4387 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4388 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4389 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4390 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4391 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4392 if (nwindow.gt.0) then
4393 allocate(winstart(nwindow)) !!el (maxres)
4394 allocate(winend(nwindow)) !!el
4395 allocate(winlen(nwindow)) !!el
4396 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4398 winlen(i)=winend(i)-winstart(i)+1
4401 if (tmax.lt.tmin) tmax=tmin
4402 if (tmax.eq.tmin) then
4406 if (nstepc.gt.0 .and. nsteph.gt.0) then
4407 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4408 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4410 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4411 ! Probabilities of different move types
4412 sumpro_type(0)=0.0D0
4413 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4414 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4415 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4416 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4417 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4418 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4419 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4421 print *,'i',i,' sumprotype',sumpro_type(i)
4422 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4423 print *,'i',i,' sumprotype',sumpro_type(i)
4426 end subroutine mcmread
4427 !-----------------------------------------------------------------------------
4428 subroutine read_minim
4431 ! implicit real*8 (a-h,o-z)
4432 ! include 'DIMENSIONS'
4433 ! include 'COMMON.MINIM'
4434 ! include 'COMMON.IOUNITS'
4435 ! character(len=80) :: ucase
4436 character(len=320) :: minimcard
4438 ! integer :: ntf,ik,iw_pdb
4439 ! real(kind=8) :: var,ene
4441 call card_concat(minimcard,.true.)
4442 call readi(minimcard,'MAXMIN',maxmin,2000)
4443 call readi(minimcard,'MAXFUN',maxfun,5000)
4444 call readi(minimcard,'MINMIN',minmin,maxmin)
4445 call readi(minimcard,'MINFUN',minfun,maxmin)
4446 call reada(minimcard,'TOLF',tolf,1.0D-2)
4447 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4448 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4449 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4450 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4451 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4452 'Options in energy minimization:'
4453 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4454 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4455 'MinMin:',MinMin,' MinFun:',MinFun,&
4456 ' TolF:',TolF,' RTolF:',RTolF
4458 end subroutine read_minim
4459 !-----------------------------------------------------------------------------
4460 subroutine openunits
4462 use MD_data, only: usampl
4465 use control_data, only:out1file
4466 use control, only: getenv_loc
4467 ! implicit real*8 (a-h,o-z)
4468 ! include 'DIMENSIONS'
4471 character(len=16) :: form,nodename
4472 integer :: nodelen,ierror,npos
4474 ! include 'COMMON.SETUP'
4475 ! include 'COMMON.IOUNITS'
4476 ! include 'COMMON.MD'
4477 ! include 'COMMON.CONTROL'
4478 integer :: lenpre,lenpot,lentmp !,ilen
4480 character(len=3) :: out1file_text !,ucase
4481 character(len=3) :: ll
4484 ! integer :: ntf,ik,iw_pdb
4485 ! real(kind=8) :: var,ene
4487 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4488 call getenv_loc("PREFIX",prefix)
4490 call getenv_loc("POT",pot)
4491 call getenv_loc("DIRTMP",tmpdir)
4492 call getenv_loc("CURDIR",curdir)
4493 call getenv_loc("OUT1FILE",out1file_text)
4494 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4495 out1file_text=ucase(out1file_text)
4496 if (out1file_text(1:1).eq."Y") then
4499 out1file=fg_rank.gt.0
4504 if (lentmp.gt.0) then
4505 write (*,'(80(1h!))')
4506 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4507 write (*,'(80(1h!))')
4508 write (*,*)"All output files will be on node /tmp directory."
4510 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4511 if (me.eq.king) then
4512 write (*,*) "The master node is ",nodename
4513 else if (fg_rank.eq.0) then
4514 write (*,*) "I am the CG slave node ",nodename
4516 write (*,*) "I am the FG slave node ",nodename
4519 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4520 lenpre = lentmp+lenpre+1
4522 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4523 ! Get the names and open the input files
4524 #if defined(WINIFL) || defined(WINPGI)
4525 open(1,file=pref_orig(:ilen(pref_orig))// &
4526 '.inp',status='old',readonly,shared)
4527 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4528 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4529 ! Get parameter filenames and open the parameter files.
4530 call getenv_loc('BONDPAR',bondname)
4531 open (ibond,file=bondname,status='old',readonly,shared)
4532 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4533 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4534 call getenv_loc('THETPAR',thetname)
4535 open (ithep,file=thetname,status='old',readonly,shared)
4536 call getenv_loc('ROTPAR',rotname)
4537 open (irotam,file=rotname,status='old',readonly,shared)
4538 call getenv_loc('TORPAR',torname)
4539 open (itorp,file=torname,status='old',readonly,shared)
4540 call getenv_loc('TORDPAR',tordname)
4541 open (itordp,file=tordname,status='old',readonly,shared)
4542 call getenv_loc('FOURIER',fouriername)
4543 open (ifourier,file=fouriername,status='old',readonly,shared)
4544 call getenv_loc('ELEPAR',elename)
4545 open (ielep,file=elename,status='old',readonly,shared)
4546 call getenv_loc('SIDEPAR',sidename)
4547 open (isidep,file=sidename,status='old',readonly,shared)
4549 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4550 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4551 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4552 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4553 call getenv_loc('TORPAR_NUCL',torname_nucl)
4554 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4555 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4556 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4557 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4558 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4561 #elif (defined CRAY) || (defined AIX)
4562 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4564 ! print *,"Processor",myrank," opened file 1"
4565 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4566 ! print *,"Processor",myrank," opened file 9"
4567 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4568 ! Get parameter filenames and open the parameter files.
4569 call getenv_loc('BONDPAR',bondname)
4570 open (ibond,file=bondname,status='old',action='read')
4571 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4572 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4574 ! print *,"Processor",myrank," opened file IBOND"
4575 call getenv_loc('THETPAR',thetname)
4576 open (ithep,file=thetname,status='old',action='read')
4577 ! print *,"Processor",myrank," opened file ITHEP"
4578 call getenv_loc('ROTPAR',rotname)
4579 open (irotam,file=rotname,status='old',action='read')
4580 ! print *,"Processor",myrank," opened file IROTAM"
4581 call getenv_loc('TORPAR',torname)
4582 open (itorp,file=torname,status='old',action='read')
4583 ! print *,"Processor",myrank," opened file ITORP"
4584 call getenv_loc('TORDPAR',tordname)
4585 open (itordp,file=tordname,status='old',action='read')
4586 ! print *,"Processor",myrank," opened file ITORDP"
4587 call getenv_loc('SCCORPAR',sccorname)
4588 open (isccor,file=sccorname,status='old',action='read')
4589 ! print *,"Processor",myrank," opened file ISCCOR"
4590 call getenv_loc('FOURIER',fouriername)
4591 open (ifourier,file=fouriername,status='old',action='read')
4592 ! print *,"Processor",myrank," opened file IFOURIER"
4593 call getenv_loc('ELEPAR',elename)
4594 open (ielep,file=elename,status='old',action='read')
4595 ! print *,"Processor",myrank," opened file IELEP"
4596 call getenv_loc('SIDEPAR',sidename)
4597 open (isidep,file=sidename,status='old',action='read')
4599 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4600 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4601 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4602 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4603 call getenv_loc('TORPAR_NUCL',torname_nucl)
4604 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4605 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4606 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4607 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4608 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4610 call getenv_loc('LIPTRANPAR',liptranname)
4611 open (iliptranpar,file=liptranname,status='old',action='read')
4612 call getenv_loc('TUBEPAR',tubename)
4613 open (itube,file=tubename,status='old',action='read')
4614 call getenv_loc('IONPAR',ionname)
4615 open (iion,file=ionname,status='old',action='read')
4617 ! print *,"Processor",myrank," opened file ISIDEP"
4618 ! print *,"Processor",myrank," opened parameter files"
4620 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4621 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4622 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4623 ! Get parameter filenames and open the parameter files.
4624 call getenv_loc('BONDPAR',bondname)
4625 open (ibond,file=bondname,status='old')
4626 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4627 open (ibond_nucl,file=bondname_nucl,status='old')
4629 call getenv_loc('THETPAR',thetname)
4630 open (ithep,file=thetname,status='old')
4631 call getenv_loc('ROTPAR',rotname)
4632 open (irotam,file=rotname,status='old')
4633 call getenv_loc('TORPAR',torname)
4634 open (itorp,file=torname,status='old')
4635 call getenv_loc('TORDPAR',tordname)
4636 open (itordp,file=tordname,status='old')
4637 call getenv_loc('SCCORPAR',sccorname)
4638 open (isccor,file=sccorname,status='old')
4639 call getenv_loc('FOURIER',fouriername)
4640 open (ifourier,file=fouriername,status='old')
4641 call getenv_loc('ELEPAR',elename)
4642 open (ielep,file=elename,status='old')
4643 call getenv_loc('SIDEPAR',sidename)
4644 open (isidep,file=sidename,status='old')
4646 open (ithep_nucl,file=thetname_nucl,status='old')
4647 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4648 open (irotam_nucl,file=rotname_nucl,status='old')
4649 call getenv_loc('TORPAR_NUCL',torname_nucl)
4650 open (itorp_nucl,file=torname_nucl,status='old')
4651 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4652 open (itordp_nucl,file=tordname_nucl,status='old')
4653 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4654 open (isidep_nucl,file=sidename_nucl,status='old')
4656 call getenv_loc('LIPTRANPAR',liptranname)
4657 open (iliptranpar,file=liptranname,status='old')
4658 call getenv_loc('TUBEPAR',tubename)
4659 open (itube,file=tubename,status='old')
4660 call getenv_loc('IONPAR',ionname)
4661 open (iion,file=ionname,status='old')
4663 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4665 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4666 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4667 ! Get parameter filenames and open the parameter files.
4668 call getenv_loc('BONDPAR',bondname)
4669 open (ibond,file=bondname,status='old',action='read')
4670 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4671 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4672 call getenv_loc('THETPAR',thetname)
4673 open (ithep,file=thetname,status='old',action='read')
4674 call getenv_loc('ROTPAR',rotname)
4675 open (irotam,file=rotname,status='old',action='read')
4676 call getenv_loc('TORPAR',torname)
4677 open (itorp,file=torname,status='old',action='read')
4678 call getenv_loc('TORDPAR',tordname)
4679 open (itordp,file=tordname,status='old',action='read')
4680 call getenv_loc('SCCORPAR',sccorname)
4681 open (isccor,file=sccorname,status='old',action='read')
4683 call getenv_loc('THETPARPDB',thetname_pdb)
4684 print *,"thetname_pdb ",thetname_pdb
4685 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4686 print *,ithep_pdb," opened"
4688 call getenv_loc('FOURIER',fouriername)
4689 open (ifourier,file=fouriername,status='old',readonly)
4690 call getenv_loc('ELEPAR',elename)
4691 open (ielep,file=elename,status='old',readonly)
4692 call getenv_loc('SIDEPAR',sidename)
4693 open (isidep,file=sidename,status='old',readonly)
4695 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4696 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4697 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4698 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4699 call getenv_loc('TORPAR_NUCL',torname_nucl)
4700 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4701 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4702 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4703 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4704 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4705 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4706 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4707 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
4708 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
4709 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
4710 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
4711 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
4712 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
4715 call getenv_loc('LIPTRANPAR',liptranname)
4716 open (iliptranpar,file=liptranname,status='old',action='read')
4717 call getenv_loc('TUBEPAR',tubename)
4718 open (itube,file=tubename,status='old',action='read')
4719 call getenv_loc('IONPAR',ionname)
4720 open (iion,file=ionname,status='old',action='read')
4723 call getenv_loc('ROTPARPDB',rotname_pdb)
4724 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4727 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4728 #if defined(WINIFL) || defined(WINPGI)
4729 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4730 #elif (defined CRAY) || (defined AIX)
4731 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4733 open (iscpp_nucl,file=scpname_nucl,status='old')
4735 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4740 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4741 ! Use -DOLDSCP to use hard-coded constants instead.
4743 call getenv_loc('SCPPAR',scpname)
4744 #if defined(WINIFL) || defined(WINPGI)
4745 open (iscpp,file=scpname,status='old',readonly,shared)
4746 #elif (defined CRAY) || (defined AIX)
4747 open (iscpp,file=scpname,status='old',action='read')
4749 open (iscpp,file=scpname,status='old')
4751 open (iscpp,file=scpname,status='old',action='read')
4754 call getenv_loc('PATTERN',patname)
4755 #if defined(WINIFL) || defined(WINPGI)
4756 open (icbase,file=patname,status='old',readonly,shared)
4757 #elif (defined CRAY) || (defined AIX)
4758 open (icbase,file=patname,status='old',action='read')
4760 open (icbase,file=patname,status='old')
4762 open (icbase,file=patname,status='old',action='read')
4765 ! Open output file only for CG processes
4766 ! print *,"Processor",myrank," fg_rank",fg_rank
4767 if (fg_rank.eq.0) then
4769 if (nodes.eq.1) then
4772 npos = dlog10(dfloat(nodes-1))+1
4774 if (npos.lt.3) npos=3
4775 write (liczba,'(i1)') npos
4776 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4778 write (liczba,form) me
4779 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4780 liczba(:ilen(liczba))
4781 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4783 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4785 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4786 liczba(:ilen(liczba))//'.mol2'
4787 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4788 liczba(:ilen(liczba))//'.stat'
4790 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4791 //liczba(:ilen(liczba))//'.stat')
4792 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4795 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4796 liczba(:ilen(liczba))//'.const'
4801 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4802 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4803 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4804 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4805 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4807 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4809 rest2name=prefix(:ilen(prefix))//'.rst'
4811 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4814 #if defined(AIX) || defined(PGI)
4815 if (me.eq.king .or. .not. out1file) &
4816 open(iout,file=outname,status='unknown')
4818 if (fg_rank.gt.0) then
4819 write (liczba,'(i3.3)') myrank/nfgtasks
4820 write (ll,'(bz,i3.3)') fg_rank
4821 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4826 open(igeom,file=intname,status='unknown',position='append')
4827 open(ipdb,file=pdbname,status='unknown')
4828 open(imol2,file=mol2name,status='unknown')
4829 open(istat,file=statname,status='unknown',position='append')
4831 !1out open(iout,file=outname,status='unknown')
4834 if (me.eq.king .or. .not.out1file) &
4835 open(iout,file=outname,status='unknown')
4837 if (fg_rank.gt.0) then
4838 write (liczba,'(i3.3)') myrank/nfgtasks
4839 write (ll,'(bz,i3.3)') fg_rank
4840 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4845 open(igeom,file=intname,status='unknown',access='append')
4846 open(ipdb,file=pdbname,status='unknown')
4847 open(imol2,file=mol2name,status='unknown')
4848 open(istat,file=statname,status='unknown',access='append')
4850 !1out open(iout,file=outname,status='unknown')
4853 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4854 csa_seed=prefix(:lenpre)//'.CSA.seed'
4855 csa_history=prefix(:lenpre)//'.CSA.history'
4856 csa_bank=prefix(:lenpre)//'.CSA.bank'
4857 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4858 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4859 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4860 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4861 csa_int=prefix(:lenpre)//'.int'
4862 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4863 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4864 csa_in=prefix(:lenpre)//'.CSA.in'
4865 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4868 write (iout,'(80(1h-))')
4869 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4870 write (iout,'(80(1h-))')
4871 write (iout,*) "Input file : ",&
4872 pref_orig(:ilen(pref_orig))//'.inp'
4873 write (iout,*) "Output file : ",&
4874 outname(:ilen(outname))
4876 write (iout,*) "Sidechain potential file : ",&
4877 sidename(:ilen(sidename))
4879 write (iout,*) "SCp potential file : ",&
4880 scpname(:ilen(scpname))
4882 write (iout,*) "Electrostatic potential file : ",&
4883 elename(:ilen(elename))
4884 write (iout,*) "Cumulant coefficient file : ",&
4885 fouriername(:ilen(fouriername))
4886 write (iout,*) "Torsional parameter file : ",&
4887 torname(:ilen(torname))
4888 write (iout,*) "Double torsional parameter file : ",&
4889 tordname(:ilen(tordname))
4890 write (iout,*) "SCCOR parameter file : ",&
4891 sccorname(:ilen(sccorname))
4892 write (iout,*) "Bond & inertia constant file : ",&
4893 bondname(:ilen(bondname))
4894 write (iout,*) "Bending parameter file : ",&
4895 thetname(:ilen(thetname))
4896 write (iout,*) "Rotamer parameter file : ",&
4897 rotname(:ilen(rotname))
4900 write (iout,*) "Thetpdb parameter file : ",&
4901 thetname_pdb(:ilen(thetname_pdb))
4904 write (iout,*) "Threading database : ",&
4905 patname(:ilen(patname))
4907 write (iout,*)" DIRTMP : ",&
4909 write (iout,'(80(1h-))')
4912 end subroutine openunits
4913 !-----------------------------------------------------------------------------
4916 use geometry_data, only: nres,dc
4918 ! implicit real*8 (a-h,o-z)
4919 ! include 'DIMENSIONS'
4920 ! include 'COMMON.CHAIN'
4921 ! include 'COMMON.IOUNITS'
4922 ! include 'COMMON.MD'
4925 ! real(kind=8) :: var,ene
4927 open(irest2,file=rest2name,status='unknown')
4928 read(irest2,*) totT,EK,potE,totE,t_bath
4931 ! AL 4/17/17: Now reading d_t(0,:) too
4933 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4936 ! AL 4/17/17: Now reading d_c(0,:) too
4938 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4941 read (irest2,*) iset
4945 end subroutine readrst
4946 !-----------------------------------------------------------------------------
4947 subroutine read_fragments
4951 use control_data, only:out1file
4954 ! implicit real*8 (a-h,o-z)
4955 ! include 'DIMENSIONS'
4959 ! include 'COMMON.SETUP'
4960 ! include 'COMMON.CHAIN'
4961 ! include 'COMMON.IOUNITS'
4962 ! include 'COMMON.MD'
4963 ! include 'COMMON.CONTROL'
4966 ! real(kind=8) :: var,ene
4968 read(inp,*) nset,nfrag,npair,nfrag_back
4970 !el from module energy
4971 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4972 if(.not.allocated(wfrag_back)) then
4973 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4974 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4976 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
4977 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
4979 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
4980 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
4983 if(me.eq.king.or..not.out1file) &
4984 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
4985 " nfrag_back",nfrag_back
4987 read(inp,*) mset(iset)
4989 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
4991 if(me.eq.king.or..not.out1file) &
4992 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
4993 ifrag(2,i,iset), qinfrag(i,iset)
4996 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
4998 if(me.eq.king.or..not.out1file) &
4999 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5000 ipair(2,i,iset), qinpair(i,iset)
5003 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5004 wfrag_back(3,i,iset),&
5005 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5006 if(me.eq.king.or..not.out1file) &
5007 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5008 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5012 end subroutine read_fragments
5013 !-----------------------------------------------------------------------------
5015 !-----------------------------------------------------------------------------
5019 ! implicit real*8 (a-h,o-z)
5020 ! include 'DIMENSIONS'
5021 ! include 'COMMON.CSA'
5022 ! include 'COMMON.BANK'
5023 ! include 'COMMON.IOUNITS'
5025 ! integer :: ntf,ik,iw_pdb
5026 ! real(kind=8) :: var,ene
5028 open(icsa_in,file=csa_in,status="old",err=100)
5029 read(icsa_in,*) nconf
5030 read(icsa_in,*) jstart,jend
5031 read(icsa_in,*) nstmax
5032 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5033 read(icsa_in,*) nran0,nran1,irr
5034 read(icsa_in,*) nseed
5035 read(icsa_in,*) ntotal,cut1,cut2
5036 read(icsa_in,*) estop
5037 read(icsa_in,*) icmax,irestart
5038 read(icsa_in,*) ntbankm,dele,difcut
5039 read(icsa_in,*) iref,rmscut,pnccut
5040 read(icsa_in,*) ndiff
5047 end subroutine csa_read
5048 !-----------------------------------------------------------------------------
5049 subroutine initial_write
5052 ! implicit real*8 (a-h,o-z)
5053 ! include 'DIMENSIONS'
5054 ! include 'COMMON.CSA'
5055 ! include 'COMMON.BANK'
5056 ! include 'COMMON.IOUNITS'
5058 ! integer :: ntf,ik,iw_pdb
5059 ! real(kind=8) :: var,ene
5061 open(icsa_seed,file=csa_seed,status="unknown")
5062 write(icsa_seed,*) "seed"
5064 #if defined(AIX) || defined(PGI)
5065 open(icsa_history,file=csa_history,status="unknown",&
5068 open(icsa_history,file=csa_history,status="unknown",&
5071 write(icsa_history,*) nconf
5072 write(icsa_history,*) jstart,jend
5073 write(icsa_history,*) nstmax
5074 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5075 write(icsa_history,*) nran0,nran1,irr
5076 write(icsa_history,*) nseed
5077 write(icsa_history,*) ntotal,cut1,cut2
5078 write(icsa_history,*) estop
5079 write(icsa_history,*) icmax,irestart
5080 write(icsa_history,*) ntbankm,dele,difcut
5081 write(icsa_history,*) iref,rmscut,pnccut
5082 write(icsa_history,*) ndiff
5084 write(icsa_history,*)
5087 open(icsa_bank1,file=csa_bank1,status="unknown")
5088 write(icsa_bank1,*) 0
5092 end subroutine initial_write
5093 !-----------------------------------------------------------------------------
5094 subroutine restart_write
5097 ! implicit real*8 (a-h,o-z)
5098 ! include 'DIMENSIONS'
5099 ! include 'COMMON.IOUNITS'
5100 ! include 'COMMON.CSA'
5101 ! include 'COMMON.BANK'
5103 ! integer :: ntf,ik,iw_pdb
5104 ! real(kind=8) :: var,ene
5106 #if defined(AIX) || defined(PGI)
5107 open(icsa_history,file=csa_history,position="append")
5109 open(icsa_history,file=csa_history,access="append")
5111 write(icsa_history,*)
5112 write(icsa_history,*) "This is restart"
5113 write(icsa_history,*)
5114 write(icsa_history,*) nconf
5115 write(icsa_history,*) jstart,jend
5116 write(icsa_history,*) nstmax
5117 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5118 write(icsa_history,*) nran0,nran1,irr
5119 write(icsa_history,*) nseed
5120 write(icsa_history,*) ntotal,cut1,cut2
5121 write(icsa_history,*) estop
5122 write(icsa_history,*) icmax,irestart
5123 write(icsa_history,*) ntbankm,dele,difcut
5124 write(icsa_history,*) iref,rmscut,pnccut
5125 write(icsa_history,*) ndiff
5126 write(icsa_history,*)
5127 write(icsa_history,*) "irestart is: ", irestart
5129 write(icsa_history,*)
5133 end subroutine restart_write
5134 !-----------------------------------------------------------------------------
5136 !-----------------------------------------------------------------------------
5137 subroutine write_pdb(npdb,titelloc,ee)
5139 ! implicit real*8 (a-h,o-z)
5140 ! include 'DIMENSIONS'
5141 ! include 'COMMON.IOUNITS'
5142 character(len=50) :: titelloc1
5143 character*(*) :: titelloc
5144 character(len=3) :: zahl
5145 character(len=5) :: liczba5
5147 integer :: npdb !,ilen
5151 ! real(kind=8) :: var,ene
5155 if (npdb.lt.1000) then
5156 call numstr(npdb,zahl)
5157 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5159 if (npdb.lt.10000) then
5160 write(liczba5,'(i1,i4)') 0,npdb
5162 write(liczba5,'(i5)') npdb
5164 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5166 call pdbout(ee,titelloc1,ipdb)
5169 end subroutine write_pdb
5170 !-----------------------------------------------------------------------------
5172 !-----------------------------------------------------------------------------
5173 subroutine write_thread_summary
5174 ! Thread the sequence through a database of known structures
5175 use control_data, only: refstr
5177 use energy_data, only: n_ene_comp
5179 ! implicit real*8 (a-h,o-z)
5180 ! include 'DIMENSIONS'
5182 use MPI_data !include 'COMMON.INFO'
5185 ! include 'COMMON.CONTROL'
5186 ! include 'COMMON.CHAIN'
5187 ! include 'COMMON.DBASE'
5188 ! include 'COMMON.INTERACT'
5189 ! include 'COMMON.VAR'
5190 ! include 'COMMON.THREAD'
5191 ! include 'COMMON.FFIELD'
5192 ! include 'COMMON.SBRIDGE'
5193 ! include 'COMMON.HEADER'
5194 ! include 'COMMON.NAMES'
5195 ! include 'COMMON.IOUNITS'
5196 ! include 'COMMON.TIME1'
5198 integer,dimension(maxthread) :: ip
5199 real(kind=8),dimension(0:n_ene) :: energia
5201 integer :: i,j,ii,jj,ipj,ik,kk,ist
5202 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5204 write (iout,'(30x,a/)') &
5205 ' *********** Summary threading statistics ************'
5206 write (iout,'(a)') 'Initial energies:'
5207 write (iout,'(a4,2x,a12,14a14,3a8)') &
5208 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5209 'RMSnat','NatCONT','NNCONT','RMS'
5210 ! Energy sort patterns
5215 enet=ener(n_ene-1,ip(i))
5218 if (ener(n_ene-1,ip(j)).lt.enet) then
5220 enet=ener(n_ene-1,ip(j))
5232 ist=nres_base(2,ii)+ipatt(2,i)
5234 energia(i)=ener0(kk,i)
5236 etot=ener0(n_ene_comp+1,i)
5237 rmsnat=ener0(n_ene_comp+2,i)
5238 rms=ener0(n_ene_comp+3,i)
5239 frac=ener0(n_ene_comp+4,i)
5240 frac_nn=ener0(n_ene_comp+5,i)
5243 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5244 i,str_nam(ii),ist+1,&
5245 (energia(print_order(kk)),kk=1,nprint_ene),&
5246 etot,rmsnat,frac,frac_nn,rms
5248 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5249 i,str_nam(ii),ist+1,&
5250 (energia(print_order(kk)),kk=1,nprint_ene),etot
5253 write (iout,'(//a)') 'Final energies:'
5254 write (iout,'(a4,2x,a12,17a14,3a8)') &
5255 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5256 'RMSnat','NatCONT','NNCONT','RMS'
5260 ist=nres_base(2,ii)+ipatt(2,i)
5262 energia(kk)=ener(kk,ik)
5264 etot=ener(n_ene_comp+1,i)
5265 rmsnat=ener(n_ene_comp+2,i)
5266 rms=ener(n_ene_comp+3,i)
5267 frac=ener(n_ene_comp+4,i)
5268 frac_nn=ener(n_ene_comp+5,i)
5269 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5270 i,str_nam(ii),ist+1,&
5271 (energia(print_order(kk)),kk=1,nprint_ene),&
5272 etot,rmsnat,frac,frac_nn,rms
5274 write (iout,'(/a/)') 'IEXAM array:'
5275 write (iout,'(i5)') nexcl
5277 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5279 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5280 'Max. time for threading step ',max_time_for_thread,&
5281 'Average time for threading step: ',ave_time_for_thread
5283 end subroutine write_thread_summary
5284 !-----------------------------------------------------------------------------
5285 subroutine write_stat_thread(ithread,ipattern,ist)
5287 use energy_data, only: n_ene_comp
5289 ! implicit real*8 (a-h,o-z)
5290 ! include "DIMENSIONS"
5291 ! include "COMMON.CONTROL"
5292 ! include "COMMON.IOUNITS"
5293 ! include "COMMON.THREAD"
5294 ! include "COMMON.FFIELD"
5295 ! include "COMMON.DBASE"
5296 ! include "COMMON.NAMES"
5297 real(kind=8),dimension(0:n_ene) :: energia
5299 integer :: ithread,ipattern,ist,i
5300 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5302 #if defined(AIX) || defined(PGI)
5303 open(istat,file=statname,position='append')
5305 open(istat,file=statname,access='append')
5308 energia(i)=ener(i,ithread)
5310 etot=ener(n_ene_comp+1,ithread)
5311 rmsnat=ener(n_ene_comp+2,ithread)
5312 rms=ener(n_ene_comp+3,ithread)
5313 frac=ener(n_ene_comp+4,ithread)
5314 frac_nn=ener(n_ene_comp+5,ithread)
5315 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5316 ithread,str_nam(ipattern),ist+1,&
5317 (energia(print_order(i)),i=1,nprint_ene),&
5318 etot,rmsnat,frac,frac_nn,rms
5321 end subroutine write_stat_thread
5322 !-----------------------------------------------------------------------------
5324 !-----------------------------------------------------------------------------
5325 end module io_config